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\^^ Abstract 

^.^ An enormous number of model chemistries are used in computational 

ri chemistry to solve or approximately solve the Schrodinger equation; each 

^ I with their own drawbacks. One key limitation is that the hardware used in 

•j computational chemistry is based on classical physics, and is often not well 

r^ i suited for simulating models in quantum physics. In this review, we focus 

C^ on applications of quantum computation to chemical physics problems. 

^ ^ We describe the algorithms that have been proposed for the electronic- 

structure problem, the simulation of chemical dynamics, thermal state 
preparation, density functional theory and adiabatic quantum simulation. 
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1 Introduction 

Controllable quantum systems provide unique opportunities for solving prob- 
lems in quantum chemistry and many-body physics that are intractable by 
classical computers. This approach is called "quantum simulation'!^ and was 
pioneered by Feynman (1982). There are two different approaches for quantum 
simulation: analog or digital. In analog quantum simulation, dedicated physi- 
cal systems are engineered to emulate the behavior of other quantum systems. 
A classic example is the use of atoms trapped in optical lattices to simulate 
the (Bose-)Hubbard model. Analog simulators are therefore special-purposed 
machines. On the other hand, digital simulation uses a universal quantum com- 
puter. Interestingly, a universal quantum computer is also capable, in principle, 
of factoring arbitrary long numbers (Shor 1997), whereas a classical computer 
is not known to be able to perform the same task. For a recent review, see e.g. 
Ladd, Jelezko, Laflamme, Nakamura, Monroe & O'Brien (2010) and Hauke, 
Cucchietti, Tagliacozzo, Deutsch & Lewenstein (2011). 

One key advantage of simulations with quantum computers over classical 
computers is the huge Hilbert space available to faithfully represent quantum 
systems. Moreover, quantum simulation avoids many problems encountered in 
classical simulation. For example, many classical algorithms relying on Monte 
Carlo methods exhibit the so-called fermion sign problem that severely damages 
the performance of the algorithm. In quantum simulation, this problem can be 
avoided by either encoding the fully anti-symmetrized wavefunction in the qubit 

^Unfortunately, the term "quantum simulation" in the community of computational physics 
refers to numerical simulation of quantum systems using classical computers. 



states, or by performing Jordan- Wigner transformations in the second quantized 
Hamiltonian and turn it into a spin Hamiltonian as first suggested by Ortiz, 
Gubernatis, Knill & Lafiamme (2001). The latter case will often result in non- 
local interaction terms, but it can still be simulated efficiently in a quantum 
computer. 

The purpose of this article is to introduce the basic concepts of digital quan- 
tum simulation and several recent developments achieved in our group. We note 
that this is by no means a comprehensive review of the whole literature in quan- 
tum simulation. We will selectively cover materials that we find most useful to 
convey an overall idea about the current status of quantum digital simulation. 
Several review articles (Zalka 1998b, Buluta & Nori 2009, Brown, Munro & 
Kendon 2010, Kassal, Whitfield, Perdomo- Ortiz, Yung & Aspuru-Guzik 2011) 
and book chapters (Kitaev, Shen & Vyalyi 2002, Stolze & Suter 2008, Nielsen 
& Chuang 2011, Williams 2010) already present a different emphasis. This re- 
view also contains some new material. Sections |2.2.2| and |2.2.3| present new 
descriptions of the simulation in the first and second quantized representations. 



respectively. Section 2.3.3 lays out a new point of view for the perturbative 
update of thermal states from smaller to bigger quantum systems, and a new 
bound for the change of a thermal state due to a perturbation. 

1.1 Quantum computational complexity and chemistry 

1.1.1 An exponential wall for many-body problems 

The theory of computational complexity studies the scaling of the resources 
necessary to solve a given problem as a function of input size. Problems are 
considered to be "easy," or efficiently solvable, if the time (or number of steps) 
for solving the problem scales as a polynomial of the input size n. For example, 
sorting a list of n items will take at most O(n^) steps. On the other hand, 
problems are considered "hard" if the scaling is exponential in n. This expo- 
nential scaling is essentially true in the worst case for almost all many-body 
problems in physics and chemistry (Pople 1999). A concise discussion of this 
point is given by Kohn (1999), where the exponential scaling of the Hilbert space 
of many-electron problems is referred to as the "Van Vleck catastrophe." The 
argument presented is as follows: if for each molecule, the accuracy to which 
one can approximate the state is (1 — e) (under a suitable metric), then for n 
non-overlapping (and non-identical) molecules, the approximation worsens ex- 
ponentially as (1 — e)". In the next subsection, we discuss the connection of 
many-body problems with computational complexity further. 

1.1.2 Computational complexity of quantum simulation 

The study of the computational complexity of problems in quantum simulation 
helps us better understand how quantum computers can surpass classical com- 
puters. It has also spurred new developments in computational complexity. For 
simplicity, computational complexity is often formulated using decision prob- 
lems. A decision problem resolves if some condition is true or false e.g. is the 



ground-state energy of the system below a certain critical value? Although the 
answer to decision problems is either "yes" or "no," one can keep asking ques- 
tions in a binary search fashion. For instance, one could attempt to determine 
in this way the ground-state energy to an arbitrarily high accuracy. 

A complexity class contains a set of computational problems that share 
some common properties about the computational resources required for solv- 
ing them. We briefly summarize below a few important examples of complexity 
classes of decision problems. 

P and NP problems 

The complexity class P contains all decision problems that are solvable in a 
polynomial time with a classical computer (more precisely, a deterministic Tur- 
ing machine). Roughly speaking, solving a problem in a polynomial time refers 
to the cases where the number of steps for solving the problem scale as a polyno- 
mial power instead of exponentially. This is considered "efficient" but, of course, 
there could be exceptions. For example, problems that scale as ©(n^*'"™) may 
take very long time to finish, compared with ones that scale exponentially as 
0(1.0001"). 

Nevertheless, from a theoretical perspective, this division allows for con- 
siderable progress to be made without considering the minutiae of the specific 
system or implementation. However from a practical standpoint, the order of 
the polynomial may be very important; especially in chemistry where an algo- 
rithm is applied to many molecules and many geometries. That said, the notion 
of polynomial versus exponential makes sense when considering Moore's "law:" 
the density of transistors in classical computers doubles every two year^ If the 
algorithm runs in exponential-time, one may be forced to wait several lifetimes 
in order for an instance to become soluble due to better classical hardware. 

Practically, common hard problems typically fall into the complexity class NP, 
which contains decision problems whose "yes" instances can be efficiently ver- 
ified to be true with a classical computer given an appropriate "solution" or 
witness. There is no doubt that P is a subclass of NP, i.e., 

PcNP . (1) 

As an example, finding the prime factors of an integer belongs to an NP problem; 
once the factors are given, then it is easy to check the answer by performing a 
multiplication. Interestingly, finding the ground state energy of the Ising model 

J2 <Ti<ji + j2<ji, (2) 

where (V, E) is a planar graph, is an NP-complete (Barahona 1982). This implies 
that if a polynomial algorithm for finding the ground state energy is found, then 



■^The exponential growth in the computational density is expected to cease sometime this 
century highlighting the importance of new methods of computation such a quantum compu- 
tation. The growth in CPU clock speed has already ceased. 



all of the problems in NP could be solved in polynomial time. In other words, 
it will imply P = NP, a result considered highly unlikely. A rigorous proof or 
disproof of this statement would constitute a significant breakthrougf(j 

It is believed, but not known with certainty, that quantum computers are not 
capable of solving all NP problems efficiently. Nevertheless, as mentioned above, 
they can solve the integer-factoring problem efficiently. It is believed that the 
complexity of integer-factoring is intermediate between P and NP (Shor 1997). 

BQP and QMA problems 

The quantum analog of P and NP problems are, respectively, the BQP (bounded- 
error quantum polynomial time) and QMA (quantum Merlin Arthur) problemtrl 
BQP is the class of (decision) problems that are solvable by a quantum com- 
puter in polynomial time. QMA is the class of (decision) problems that can 
be verified by a quantum computer in polynomial time. Like NP-problems, 
the QMA class covers many problems that are important to physics and chem- 
istry (Liu, Christandl & Verstraete 2007b, Schuch & Verstraete 20G9a, Wei, 
Mosca & Nayak 2010). For example, the ground-state problem of Hamilto- 
nians involving local interaction terms is known to be QMA-complete (Kitaev 
et al. 2002, Kempe, Kitaev & Regev 2006a). For more discussion on topics of 
computational complexity and quantum simulation, readers may find the fol- 
lowing references useful: Aharonov & Naveh (2002), Rassolov & Garashchuk 
(2008), Aaronson (2009), Osborne (2011) and Aaronson (2011). 

The key point here is that so far it is not known whether quantum computers 
can solve NP and QMA problems efficiently. In fact, many attempts (see e.g. 
Young, Knysh & Smelyanskiy (2008), Poufin & Wocjan (2009a), Poulin & Woc- 
jan (2009b), Young & Smelyanskiy (2010), Bilgin & Boixo (2010), Yung, Nagaj, 
Whitfield & Aspuru-Guzik (2010), and Hen & Young (2011)) show that expo- 
nential resources are required to solve problems in these classes. Nevertheless, 
many problems in physics and chemistry do exhibit symmetries and structures 
that we could exploit to construct efficient quantum simulation algorithms. This 
is the main theme of the discussion in the rest of the paper. 

1.2 Basic quantum algorithms for digital quantum simu- 
lation 

Digital quantum simulation cannot be easily understood without a detour into 
the basics of quantum algorithms. Quantum algorithms are procedures for ap- 
plying elementary quantum logic gates to complete certain unitary transforma- 
tions of the input state. The quantum computer state is usually written in 
terms of qubits (two-level systems) . In the two-dimensional Hilbert space of a 



•^P vs NP is one of the Millennium Problems of the Clay Mathematics Institute 
[http : //www. claymath.org/niillenniuni/P_vs_NP/] 

*More precisely, BQP is analogous to the classical complexity class BPP, which refers to 
problems that can be solved with randomized algorithms in a classical computer in polynomial 
time, subject to a bounded error probability. 



single qubit, we label the upper and lower eigenstates of a^ as |0) and |1). Note 
that the choice of u^ as the computational basis is arbitrary. This is called the 
computational basis, and the matrix representation of operators and states are 
written in this basis unless otherwise stated. The unitary transformations of 
the qubits may be visualized using quantum circuit diagrams introduced later 
to explain some of the more complex quantum algorithms. 

It is known that any unitary gate can be decomposed into some sets of univer- 
sal quantum logic gates that contains single- and two-qubit operations (Nielsen 
& Chuang 2011). The first gate of interest is the single-qubit Hadamard trans- 
formation defined (in the computational basis) as 



\/2 



1 1 
1 -1 



The Hadamard gate transforms between the a'' basis and the cr^ basis (|±) = 
(|0)±|1))/a/2) and will be used throughout the article. A second gate of interest 
is the CNOT (controlled not) gate which is a non-trivial two-qubit gate. It 
leaves one input qubit unchanged and acts with g'-^ = |0)(1| + |1)(0| on the 
second qubit when the first qubit is in the state |1). The first qubit is called 
the control and the NOT operation is applied coherently when the control qubit 
is in a superposition of computational basis states. Symbolically, the gate is 
written as CNOT= |1)(1| (g) cr^ -h |0)(0| (g) /. The Hadamard and CNOT gates 
are not universal for quantum computation, and in fact quantum algorithms 
with only these gates can be simulated efficiently classically as shown by the 
Knill-Gottesman theorem (Nielsen & Chuang 2011). Therefore, this gate set 
must be augmented by other single qubit gates which can always be expressed 
by single-qubit rotations, Rx, Ry, and Rz where Rx is defined as exp[— icr'^6'/2] 
for real angle 9. 

There are two elementary algorithms, namely quantum Fourier transform 
(QFT) and phase estimation algorithm (PEA), that play important roles in 
many applications in quantum simulation. We turn our attention to them now. 

1.2.1 Quantum Fourier transform (QFT) 

Given a vector with N elements {xq^xi, ..,xn-i), in classical computation, the 
discrete Fourier transform outputs another vector of iV numbers {yo, yi, ■-, 2/JV-i) 
through the following relation: 

In quantum computation, for any given quantum state, \(f>) ~ '^x=o 'Pi^) \^)j 
the goal of the quantum Fourier transform C/qft is to perform the following 
unitary transformation: 

N-l 
fc=0 



where (f> (fc) — (1/\/N) X^x Jo '^ (a:)e^'^*^'^'^ is the Fourier-transform of the func- 
tion (t){x) (compare with Eq. (IS])). Due to the hnearity of C/qft, it is sufficient 
to consider the transformation of the basis vectors such that 

N-l 

t/QFT|x) = ^^e2-'=/^|A:) . (5) 

For a system of n qubits, the number of gates required for such a transforma- 
tion is 0{'n?) (Nielsen & Chuang 2011). For the classical case (see Eq. ([3|), one 
will require 0(n2") gates to complete the same transformation, e.g. with Fast 
Fourier transform (FFT). This may seem to suggest that quantum computers 
are exponentially more efficient in performing the task of discrete Fourier trans- 
formation. However, the caveat is that one cannot directly compare QFT with 
the classical FFT. The reason is that if we want to obtain a particular Fourier- 
transform coefficient, say (j) (k), from the quantum state in Eq. W, it would still 
require exponentially many steps to extract the information (phase and ampli- 
tude), e.g. through quantum state tomography where many measurements are 
used to analyze the state (Nielsen & Chuang 2011). 

Nevertheless, QFT is essential in many applications in digital quantum sim- 
ulation. As we shall see in section [2. 2. 2[ it allows us to simulate the time dy- 
namics of particles efficiently by moving between the position and momentum 
representations. Another important application of the QFT is phase estimation, 
discussed next. 

1.2.2 Phase estimation algorithm (PEA) 

The phase estimation algorithm C/pEA is an essential component for many quan- 
tum algorithms for quantum simulation, as well as the celebrated factoring 
algorithm (Shor 1997). Loosely speaking, the PEA can be considered as a real- 
ization of the von Neumann measurement scheme (without the last projective 
measurement) in the eigenvalue basis \ak) of any Hermitian observable A (e.g. 
Hamiltonian H). More precisely, if we prepare a register of m ancilla qubits 
initialized in the state |000...0), then for any given state \(f>) = X^/c ^fe \o-k), we 
have 

C/pEA \^) I000...0) « ^ Cfc Iflfe) \Ak) , (6) 

k 

where, for the moment, we assume that the ^^'s are the tti- integer-digit repre- 
sentation (i.e., Ak £ {0, 1, 2, .., 2™— 1}) of the eigenvalue of A. xThe projective 
measurement cannot be implemented perfectly in general (hence the ~ sym- 
bol). We will see where the errors come from as we go through the details of 
the algorithm below. 

Suppose that we are given an eigenstate \ak) of the Hermitian observable 
A. The goal of PEA is to determine Ak, given that we are able to simulate a 
unitary operator W where 

W^|a,)=e2"^M«fe) , (7) 



and (pk = ^fc/2™. The first step of the PEA is to apply Hadamard gates to each 
of the ancilla qubits. This results in an equal superposition of states 

|S> = ^|'W (8) 

where a: is a ?n-digit binary number. Then, taking each ancilla qubit j as a 
control qubit, we apply the controlled- W^^~^ gate to the state IS*) \ak); this 
effectively performs the following operation: 

\x)\ak)^\x)W^\ak) . (9) 

Of course, from Eq. (I?]), the right-hand side gives only a phase factor, namely 
exp {2Trix(f>k)- The resulting state is 



E 



e2--^'= |.t) \ak) . (10) 

\ ^ " x=0 / 

Comparing this state with that in Eq. Bl , and assuming the special cases where 
the phase angle <j)x can be expressed exactly by m binary digits, the application 
of the inverse of the quantum Fourier transform C/qft will convert the state in 



Eq. ( 10 1 into the following state 

\Ak)\aK) . (11) 

Since the unitary operator C/pEA is linear, the procedure applies to any initial 
state. For this particular case, where A^s are integers, we have shown that 
PEA is effectively a projective measurement as advertised in Eq. ^. 

For the general case, where the A}^^s are real numbers, the corresponding 
(jife's will have precision beyond 1/2™; this is the source of the errors in the 
expression of Eq. ([6| . The overall error decreases when we increase the number of 
ancilla qubits and perform several QFTs in parallel (we refer to Kaye, Laflamme 
& Mosca (2007) for a detailed error analysis). More precisely, if we want to 
achieve a p-hii precision of 0^, with an error less than e, one will need more than 
m = p-l-log (2 -t- l/2e) ancilla qubits. In general, implementing the operator W^ 
requires k times as many resources as those needed for simulating W . Therefore, 
the scaling of the quantum gates of PEA grows exponentially when we increase 
the precision p of the phase measurement. This result is consistent with that of 
the general sampling theory in classical signal processing, where the precision 
of the Fourier spectrum 5uj goes as the inverse of the total time T sampled, i.e., 
6ui ~ 0{\/T). This is because the cost of the quantum simulation is proportional 
to T, and T grows exponentially with the number of bits of precision. 

2 Digital quantum simulation 

2.1 Overview 

Broadly speaking, the steps involved in carrying out a digital quantum simu- 
lation consist of three parts: state preparation, time evolution, and measure- 



ment of observables. Measurement of Hermitian observables can be achieved via 
the phase estimation method (Abrams & Lloyd 1999, Jaksch & Papageorgiou 
2003, Knih, Ortiz & Somma 2007) described before. Other apphcations (Lidar 
& Wang 1999, Wu, Byrd & Lidar 2002, Somma, Ortiz, Gubernatis, Knill & 
Lafiamme 2002, Somma, Ortiz, KniU & Gubernatis 2003, Byrnes & Yamamoto 
2006, Kassal & Aspuru-Guzik 2009) or quantities of physical interest such as 
the partition function (Master, Yamaguchi & Yamamoto 2003, Wocjan, Chiang, 
Nagaj & Abeyesinghe 2009) , can be obtained through variants of the methods 
employed in state preparation and time evolution, and we will skip them in this 
review. Below we give an overview of state preparation and simulation of time 
evolution. It turns out that many methods of state preparation also depend 
on the time evolution itself. Therefore, we will first cover the methods of time 
evolution before state preparation. 

2.2 Simulation of time evolution 

The simulation of the time evolution of quantum state \^p) under Hamiltonian 
H according to the Schrodinger's equation {h = 1), 

z||^)=i/(i)|V) , (12) 

is one of the key applications of quantum computation. If, for example, the 
time-evolution operator 

U{t) = exp{~iHt) (13) 

can be simulated efficiently, then the eigenvalues of H can be obtained through 
the phase estimation algorithrrr] As mentioned in the introduction, Feynman 
(1982) investigated the possibility of simulating quantum systems using another 
quantum system, and conjectured that there existed a class of universal quan- 
tum simulators that evolved under a Hamiltonian with local interactions. This 
conjecture was justified by Lloyd (1996), who argued that any Hamiltonian 

771 

H = Y,H, (14) 

i=l 

which can be decomposed into m local terms {Hi} can be simulated efficiently 
by a universal quantum computer. Each Hi term acts on at most k qubits 
(or quantum subsystems). The key idea is based on the Trotter splitting or 
"trottcrization" of all non-commuting operators. 



^~^Ht 



f~iHit/n^-iH2t/n^^^^-iH^t/n\ /-^^\ 



where the approximation can be made arbitrarily tight by refining the time- 
slicing, i.e., increasing n. 



^Moreover, it can also be exploited for quantum cooling (see section 2.4 i 



There exist higher order approximations (Suzuki- Trotter formulas) which 
reduce the error even further. For instance, the second-order approximation is 
given by 



+ 0(t(At)2) (16) 

A quantum circuit on n qubits which approximates U{t), with error at most e, 
is efficient if the number of one- and two-qubit gates involved is polynomial in 
the scaling of the problem, i.e., poly(n,T, 1/e) with r = t/||-ff||. 

2.2.1 Suzuki- Trotter formulas 

We now briefly review the use of Suzuki- Trotter formulas in quantum simu- 
lation for time-independent sparse Hamiltonians, providing an introduction to 
the quantum simulation literature. Continuing the work of Lloyd (1996), works 
by Aharonov & Ta-Shma (2003) and Childs (2004) show that black-box sparse 
Hamiltonians are too efficiently simulatable. Sparsity here means that the num- 
ber of elements per row is bounded by some polynomial of n, while the dimension 
of the Hilbert space is Z? = 2". It is also required that each matrix element 
can be retrieved efficiently. Ref. (Aharonov & Ta-Shma 2003) used a color- 
ing scheme to decompose the Hamiltonian into a sum of 2 x 2 block diagonal 
matrices. This coloring scheme has been updated in several references (Berry, 
Ahokas, Cleve & Sanders 2006, Berry & Childs 2009, Childs & Kothari 2011). 
The coloring scheme and blackbox simulation will not be discussed further. 

Berry et al. (2006) were the first to approach the general problem of simulat- 
ing non-commuting Hamiltonians by using higher order Suzuki- Trotter formulas. 
Papageorgiou & Zhang (2010) returned to this issue and their contributions will 
be discussed later. The important results of Berry et al. (2006) are 

1. the use of higher order Trotter-Suzuki decompositions to bound the num- 
ber of non-commuting exponentials, N^xp, necessary to carry out a simu- 
lation for some amount of time t, 

2. a proof of a no-go theorem for sub-linear black-box simulation and 

3. improvements upon the coloring scheme of (Aharonov & Ta-Shma 2003) 
and (Childs 2004) for black box simulation of sparse Hamiltonians. 

The simulations in this review are concerned with the first two results and they 
will be explained in more detail after describing the Suzuki- Trotter formulas. 

M. Suzuki has studied and extended the Trotter formula essentially con- 
tinuously since 1990 and this work was reviewed in Hatano & Suzuki (2005). 
The recursive formulas introduced by Suzuki define a fractal pattern where a 
combination of forward and backward propagation leads to an improved approx- 
imation of the desired exponential. Suzuki defines higher order Trotter formulas 



10 



in a recursive way. Beginning with the spht operator formula, g^^/^e^^e"*^/^, 
for m operators, the foUowing scries of equations were derived: 



52(x) = n^'^'' l[e'^^^\ (17) 

\fe=l / \fe=m / 

S4{x) = S2{z2xfS2{{l-4z2)x)S2{z2xf (18) 



S2kix) = Sk{zkxfSk{{l-4,Zk)x)Sk{zkx)'^ (19) 

The values of the constants {zj} are selected so that 5*2^ is correct through 2j*^ 
order and it can be shown (Hatano & Suzuki 2005) that Zi = (4 - 4^/'^*^^))^^ 
If there are m non-commuting Hamiltonians, then the first order approximation 
takes TO = -/Veaip and for the split operator formula, 5*2 , the number of exponen- 
tials is 2m — 1. In general, 2(to — 1)5*^"^ -I- 1 exponentials are used for the S2k 
approximant. 

For the fc*'* order Suzuki- Trotter, with m Hamiltonians in the sum, and error 
tolerance given by e. Berry et al. (2006) gives a bound on the number of expo- 
nentials by bounding each order of the Suzuki- Trotter formula. Papageorgiou 
& Zhang (2010) presented an improvement by noting that the relative ratio of 
Hamiltonian norms is also important. The main idea is that if some of the 
Hamiltonians have very small weight, then their exponentials can be effectively 
ignored. 

The optimal order of Trotter decomposition, fc*, is determined by selecting 
the best compromise between time-step length and a decomposition using more 
exponentials. In Berry et al. (2006) this was worked out for unstructured sparse 
Hamiltonians N^xp > ||-ff||i. The lower bound on the generic cost of simulating 
an evolution was by contradiction, and relied on the lower bounds to quantum 
mechanical problems based on the polynomial method (Beals, Buhrman, Cleve, 
Mosca & de Wolf 1998). This bound could be violated given sub-linear simu- 
lation time. In a departure from the methods discussed so far, Childs (2010) 
used a quantum walk based approach to push the scaling closer to linear in 
the reweighed time and Raesi, Wiebe & Sanders (2011) looked at designing 
quantum circuits for quantum simulation. 

For problems in chemistry, it is more natural to represent Hamiltonians in 
terms of first- and second-quantized forms. In the following, we will describe 
how to exploit the special structure of molecular Hamiltonians to simulate the 
time dynamics in quantum computers. 

2.2.2 First-quantized representation 

In the first-quantized form, the non-relativistic molecular Hamiltonians H de- 
composes in a kinetic T and potential V terms, i.e., 

H = T + V . (20) 
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The kinetic term includes the contribution from the nuclei and electrons sepa- 
rately, 

^-i:^v.'-i:|;v? , (21) 

where Mi is the mass of the nucleus «, and nie is the electron mass. The potential 
energy term comprises of the Coulomb interaction among the nuclei, among the 
electrons, and between the nuclei and electrons. Explicitly: 

V'(r,R) = -^V-M^ + -^V,-^--^V-^^, (22) 



e2 y^ Z,Zj 1 e' y- 1 
Aneo ^ |R, - Rj 1 Atteo ^ jr^ - r^ | 


e^ y^ Z, 



where e is the electric charge, and Zi is the charge of nuclei i. The coordinates 
of nuclei i and electron j are denoted by R^ and Vj.Wc will use the notation 
r = (ri, r2, r3...) (and similarly for R). We also ignore the spin degrees of 
freedom, which can be incorporated easily. 

The general wavefunction can be represented in the position basis as 

|*)=;^*(r,R)|rir2r3...)|RiR2R3...> , (23) 

r,R 

where each electronic or nuclear coordinate is represented on its own grid over 
m qubits resulting in a total of Bm qubits to represent the state of B particles. 
Note that the grid encoded in m qubits has 2™ points. The complex wave- 
function \1/ (r, R) in addition to being properly normalized must also be anti- 
symmetrized (or symmetrized for Bosons). Abrams & Lloyd (1997) and Ward, 
Kassal & Aspuru-Guzik (2009) consider the necessary anti-symmetrization pro- 
cess for fermions in first quantization. 

To simulate the dynamics (Zalka 1998b, Wiesner 1996, Kassal, Jordan, Love, 
Mohseni & Aspuru-Guzik 2008), we note that although the kinetic and potential 
terms do not commute with each other, both can be represented as diagonal 
operators in momentum and position basis respectively. By using the quantum 
Fourier transform J/qft, it is natural to decompose the time evolution as 



^-^Ht 



('^QFTe"^"/"C^QFTe-^^*/")" . (24) 



In fact, this method is known as the split-operator method (Feit, Fleck Jr & 
Steiger 1982, Kosloff 1988). Higher-order Suzuki- Trotter formulas can also be 
applied, as described before. This method was applied to quantum computing in 
a number of works (Wiesner 1996, Zalka 1998a, Zalka 1998b, Strini 2002, Bcnenti 
& Strini 2008, Kassal et al. 2008). 

In the context of quantum computing, it remains to find a method to induce 
a coordinate-dependent phase factor such that 

|rir2r3...)|RiR2R3...) ^ e-'^f""'^)^* |rir2r3...) IR1R2R3...) , (25) 



12 



where 5t = t/n, and similarly for the kinetic term in the Fourier basis. An 
efficient methocjjis implicitly described in the book Kitaev et al. (2002) (pages 
131-135), which was further developed and adapted to the chemistry context 
by Kassal et al. (2008). We sketch the idea here for completeness. First, we 
will assume that the potential energy term is rescaled to become dimensionless, 
and projected into a range of integer values such that < V^ (r,R) < 2™ — 1, 
where m should be sufficiently large to allow appropriate resolution of V (r, R) 
in the integer representation. Next, we define a more compact notation |r, R) = 
|rir2r3...) IR1R2R3...), and an algorithmic operation A to be performed in the 
position basis: 

y^|r,R)|s) ^ |r,R)|s®F(r,R)) , (26) 

where |s), s ~ 1,2,3,..., is a quantum state of m ancilla qubits, and © is 
addition modulo 2™. Suppose now that the ancilla qubits are initialized in the 
following state: 

|q)^-L^e2-/M|,^ (27) 



s=0 

where M = 2™. This state is the Fourier transform of ll). Then the desired 



phase generating operation described in Eq. (25) can be achieved using con- 
trolled a^ rotations after applying A to the state |r, R) |q). A similar procedure 
is applied to the kinetic term to complete the Trotter cycle. 

An alternative approach to implement the controUed-phase operation de- 



scribed in Eq. ( 25 1 is the following: first include a register of qubits initialized 



as |0). Then in a similar (but not identical) way as that described in Eq. (26), 
we define the operation 

^|r,R)|0) ^ |r,R)|y(r,R)5t) , (28) 

where we used ®V {r,'R) St ^ V (r, R) St. The state \V (r, R) St) is the binary 
representation {xiX2X3...Xm} defined through the following equality, 

m 

V {r,'R) St = 2t: x0.xiX2X3...x„r = 2TT^^ . (29) 

fc=i 

Now, we can decompose the overall phase as follows, 

^-iV{r,R)St ^ ^-i2Trxi/2^-t27TX2/2^ g-i27rK™/2" /^qX 

This decomposition can be achieved through the application of m local phase 
gates Rk = |0) (0| -I- exp (— 27rz/2'^) |1) (1| for each ancilla qubit. This approach 
requires the ancilla to be un-computed (i.e. the inverse of the operation in 



Eq. (28)) in the last step. 



^An alternative method was proposed by Benenti & Strini (2008), but it scales exponen- 
tially with the number qubits. 
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2.2.3 Second-quantized representation 

The first-quantization method is universahy apphcable to any molecule. The 
shortcoming is that it does not take into account the physical symmetrization 
properties of the underlying quantum system. When a suitable set of basis 
functions is employed, the size of the problem can be significantly reduced. This 
is known as the second-quantization approach in quantum chemistry, which can 
be extended for quantum simulation. 

Most studies on quantum simulation based on first quantization methods use 
grids to represent wave functions, while works employing second quantization 
methods generally use atomic or molecular orbitals as a basis set for the wave 
functions. We will take the later approach here. Nevertheless, the choice of basis 
is not the key difference between the first and second quantization. Indeed, a 
basis set of delta functions (or approximations to delta functions) could be used 
to represent a grid within second quantization. On the other hand, the storage 
of the same wave function is very different in second and first quantization. For 
example, a two-particle wave function with the first particle at site i and the 
second at site j, is represented as | coords )| coord j) in first quantization, and as 
\0 ■ ■ ■ li ■ ■ ■ Ij ■ ■ ■ 00) in second quantization. 

The starting point of the second-quantization approach (Aspuru-Guzik, Du- 
toi. Love & Head-Gordon 2005, Wang, Kais, Aspuru-Guzik & Hoffmann 2008, 
Whitfield, Biamonte & Aspuru-Guzik 2011) is the Born-Oppenheimer approxi- 
mation, where the nuclear coordinates R are taken to be classical variables. This 
allows us to focus on the electronic structure problem. Ignoring the nuclear ki- 
netic and the nuclear-nuclear interaction terms, the molecular Hamiltonian in 
Eq. ( 20 1 can be expressed as 

if = ^ hpgalaq + ^^^ hpqrsalalaras . (31) 

pq pqrs 

where the fermionic creation operator aj, creates an electron in the p mode from 
the vacuum, i.e., aj, |vac) = \py Denote Xp(r) as the single- particle wavefunc- 
tion corresponding to mode prj Then, the explicit form for the single-electron 
integrals is given by 

,.„.-/*X;w(|-V' + 3|;I:jS^)x,W . (32) 

and the electron-electron Coulomb interaction term is, 

e^ f x;(ri)x;(r2)Xr(r2)x. (i-i) ,^^, 

hpqrs = -. / aridr2-^^ ■ . (33) 

47reo J |ri-r2| 

These integrals have to be pre-calculated with classical computers before encod- 
ing them into the quantum algorithms. If we keep k single-particle orbitals, then 



'^Here r refers to the coordinates of one particular electron. 
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there are 0{k^) terms. More details of the formaUsm of second-quantized elec- 
tronic structure theory in the Born-Oppenhcimcr approximation can be found 
in Helgaker, Jorgensen & Olsen (2000). 

To simulate time dynamics in a quantum computer, we can apply the same 



Trotterization idea described above (see Eq. (15)), and simulate separately the 
terms 

exp{—ihpqa'paqSt) and exp{—ihpqrsapa'garasSt) . (34) 

Since the simulation of every single exponential term in a quantum computer is 
costly, due to error-correction overheads as discussed in Clark, Metodi, Gasster 
& Brown (2009), one simplification we can make is to group the terms of single- 
particle terms into two-particle terms. This is possible for electronic problems 
with a fixed number N of electrons. Consider any A^-fermionic state, then the 
identity operator I^ is equivalent to a summation of the following single-body 
number operators, 

{l/N)Y,4as ^ In , (35) 

which means that we can write 

alaq = ^^— Y Yl ap^^a^Og . (36) 

s 

The last equation is a sum of two-electron terms, and can be absorbed into 
the pre-computed values of hpgrs- Now, denoting the new values as hpqrs, the 
Hamiltonian H reduces to 

F = - 2J hpqrsO'pa^aras ■ (37) 

pqrs 



Therefore, we are left only with simulating the two-body term in Eq. (34). 

One challenge we need to overcome is the fermionic nature of the operators 
al and a^, which comes from the anti-symmetrization requirement of fermionic 
wavefunctions. A first step to overcome this challenge is to map the occupation 
representation to the qubit configuration. Explicitly, for each fermionic mode 
J, we represent the qubit state |0) ■ = \\.) , as an unoccupied state, and similarly 
|1) = It) as an occupied state. To enforce the exchange symmetry, we apply 
the Jordan- Wigner transformation (Ortiz et al. 2001, Whitfield et al. 2011): 

n <^ <rj and a, = I n < < > (38) 



4 



^m<j I \^J^<J 



where 



a± = (cr^ ± iay)l2 . (39) 

By using Eq. ( 38 ) and ( 39 ) , we can now write the fermionic Hamiltonian in 



Eq. (37) as a spin Hamiltonian involving products of Pauli matrices {o^ , cr'^ ,a^}: 



Jospin — / ^ / ^ 9pqrs"pqrs<^p'^q'^r^ s i V^^j 

pqrs abed 
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where the set of indices {p, q, r, s} is summed over the ferraionic modes, and 
{a,b,c,d} is either x or y. The operator Opqrs keeps track of the ctz's; for 
example, if p > q > r > s, we then have 

n -n X ( n -I I ■ (41) 

The punchhne here is that the Hamiltonian becomes a polynomial sum of prod- 
ucts of spin operators, and each operator is locally equivalent to tr^. Therefore, 
the non-trivial part of simulating the time dynamics of the fermionic Hamilto- 
nian is to simulate the non-local interaction terms of the following form: 

exp(-i5crW^...CT^(5i) , (42) 

where g is some constant. This can be achieved by a series of controUed-NOT 
together with a local operation (see e.g. Figure 4.19 of Nielsen & Chuang 
(2011)), or the phase generating method similar to the one described in the 



previous section (cf. Eq. (25)). The explicit circuits for simulating the time 



evolution operators can be found in Whitfield et al. (2011). 

2.2.4 Open-system dynamics 

In quantum mechanics, the time evolution dynamics of a closed system is al- 
ways described by a unitary transformation of states, U (t) pW (t). However, 
non-unitary dynamics occurs when the dynamics of the system of interest S is 
coupled to the environment B, as in, 

psit)=TTB[Uit)psBUHt)] . (43) 

After some approximations this evolution can often be described by a (Marko- 
vian) quantum master equation in Lindblad form (Breuer & Petruccione 2002, 
Lindblad 1975, Gorini, Kossakokowski & Sudarshan 1976), 



di' 



p^(t) = -i[iJ,,p,]-|-^m„^( A„p,,A)j -t- A„,p,A^ ) , (44) 



Q,/3 



where H^ is the system Hamiltonian, ma^ is a positive matrix, and Aq, is a linear 
basis of traceless operators. This quantum master equation is relevant in many 
physical, chemical, and biological processes at finite temperature (Mohseni, 
Rebentrost, Lloyd & Aspuru-Guzik 2008, Rebentrost, Mohseni, Kassal, Lloyd & 
Aspuru-Guzik 2009). Further, this equation has many applications in quantum 
information processing, including preparing entangled states (from arbitrary 
initial states) (Kraus, Biichler, Diehl, Kantian, Micheli & ZoUer 2008, Krauter, 
Muschik, Jensen, Wasilewski, Petersen, Cirac & Polzik 2011, Muschik, Polzik 
k Chac 2011, Cho, Bose & Kim 2011, Miiller, Hammerer, Zhou, Roos & ZoUer 
2011), quantum memories (Pastawski, Clemente & Cirac 2011), and dissipative 
quantum computation (Verstraete, Wolf & Ignacio Cirac 2009). It has been 
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shown that the quantum master equation can be simulated by a unitary quan- 
tum circuit with polynomial resource scaling (Bacon, Childs, Chuang, Kempe, 
Leung & Zhou 2001, Khesch, Barthel, Gogolin, Kastoryano & Eisert 2011). The 



basic idea is as follows: we first re- write the master equation (Eq. (44)) in the 
form, 

j^pAt)=C{ps) , (45) 

where £ is a super-operator. Similar to the unitary dynamics, we can define the 
super-operator version of the propagator /C (ii,to) through the relation, 

Ps{ti)=IC{h,to){ps{to)) (46) 

for all values of time ti > t^. Suppose we consider a finite time interval T, 
which can be divided into m small time intervals Ai, i.e., T — mAt. Then 
similar arguments (Kliesch et al. 2011) based on Trotterization show that the 
following approximation, 

/C(r)«/C(At)/C(At)/C(At).../C(Ai) , (47) 



indeed converges when the division size goes to zero, i.e.. At — >■ 0. The remaining 
part of the argument is to show that each of the small-time propagator terms 
/C(Ai) can be simulated efficiently with a quantum circuit. This is generally 
true if the superoperator £ is a finite (polynomial) sum of local terms (Bacon 
et al. 2001). 

2.3 State preparation 

We have discussed how quantum dynamics can be simulated efficiently with 
a quantum computer, but we have not yet discussed how quantum states of 
physical or chemical interest can be initialized on the quantum computer. In 
fact, both thermal and ground states of physical Hamiltonians can be prepared 
by incorporating the methods of simulating the time dynamics, as we shall 
explain later in this section. 

We first consider a strategy to prepare quantum states that can be efficiently 
described by some integrable general function, e.g., a Gaussian wave packet. 
Before we provide a general description, it may be instructive to consider the 
case of creating a general (normalized) two-qubit state, 

/oo |00) + /oi |01) + /lo |10) + /n 111) (48) 

from the initial state |00). First of all, we will assume that all the coefficients 
/ij's are real numbers, as the phases can be generated by the method described 



in Eq. ( 25 ) . Now, we can write the state in Eq. ( 48 1 as 



<?o|0)^(|f|0) + ^|l))+gi|l)®(ff |0) + ^|1)) , (49) 
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a//qo + /oi is the probability to find the f irst qubit in the state |0), 



\/fiQ + fii ■ The form in Eq. ( 49 1 suggests that we can 



where go 

and similarly for gi 

use the following method to generate the general state of Eq. (48) from |00). 

-^ go |0) + gi\l), to the first qubit. The 



1. Apply a rotation, such that |0 
resulting state becomes, 

(5o|0)+gi|l))|0) 

2. Perform the following controlled operation: 



|0)-^ 






|0) 



/x] 

Sx 



(50) 



(51) 



where 



{0,1}. 



The final state is exactly the same as that in Eq. ( 49 1 or Eq. ( 48 1 



Consider, more generally, the preparation of the following n-qubit quantum 
state (Zalka 1998b, Grover & Rudolph 2002, Kaye & Mosca 2004, Ward et al. 
2009): 



2"-l 

E 

x=0 



/(^)k) 



(52) 



Here again we will assume that f{x) is real. We can image that this is the 
wavefunction of a particle in ID. The first qubit describes whether the particle is 
located in the left half |0) or right half |1) of the line divided by L = 2" divisions. 
The first step is therefore to rotate the first qubit as cos^o |0) + sin 611 |1), where 



cos dn = 



0= E 

0<x<L/2 



fi^Y 



(53) 



represents the probability of locating the particle at the left side, i.e. < a; < 
L/2. The next step is to apply the following controlled rotation: 



|0) 



)(S^lo) 



COS 6x1 

COS 9x 



ii: 



(54) 



where 



COS^^oo 



a<x<L/4 L/4<x<L/2 



represents the probability for finding the particle in the '00' division (0 < < L/A) 
and the '01' division {L/A < < L/2) respectively; an analogous arguments ap- 
ply for the '10' and '11' divisions. In the remaining steps, similar controlled 



operations described in Eq. (54 1 are applied, which depend on the division of 



the controlling qubits. The 9 rotation angles have to be calculated explicitly. It 
is therefore necessary that the function f{x) is efficiently integrable (Grover & 
Rudolph 2002). This is expected, as otherwise such a simple algorithm would be 
able to solve the random-field Ising spin models and other NP-complete prob- 
lems. We will cover the creation of thermal states later. In the next section, we 
will consider methods for preparing ground states. 
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Figure 1: Example for the state preparation method. The space is divided in 
L = 8 divisions, the '0' division refers to the left half of the space, (0 < x < L/2), 
and similarly for the '1' division. Finer resolution is achieved by increasing the 
number of labeling digits. 

2.3.1 Preparing ground states 

Phase-estimation based methods 

Finding ground states of classical Hamiltonians, e.g. a random-field Ising model, 
is known to be NP-hard. Therefore, it is not expected that a quantum com- 
puter would be able to solve it efficiently in general. Furthermore, preparing the 
ground-state of a general quantum Hamiltonian H is even more challenging as 
both eigenvalues and eigenvectors are required to be obtained, and this problem 
belongs to the QMA complexity class, the quantum analog of NP. Fortunately, 
many problems in physics and chemistry exhibit structures and symmetries that 
allow us to arrive at solutions that are approximation of the exact solutions; for 
example, the BCS wavcfunction related to superconductivity and superfluidity, 
and the Laughlin wavcfunction related to the fractional quantum Hall effect 
(FQHE) , both provide good predictions for the corresponding many-body prob- 
lems. The quality of other approximated solutions, such as the mean-field or 
Hartree-Fock approximation, may vary from problem to problem. 

The quality of the approximated solution (or trial solution) |V't) can be 
quantified by the fidelity F defined by 



F=\{eoHT)\' 



(56) 



where |eo) is the target ground state (assumed unique) of the Hamiltonian H of 
interest. The physical meaning of F is that if one can implement a projective 
measurement {\ek) {ek\} in the eigenvector basis {|efe)} of H to the trial state 
lipr), then the probability of getting the ground state |eo) is exactly equal to 
F, and can be implemented with the phase estimation algorithm (Abrams & 
Lloyd 1999). A similar procedure can produce low energy eigenstates even if 
there is no gap (Poulin & Wocjan 2009a). 

With the methods of the previous paragraph, if the projection on the ground 
state fails, the initial approximation must be reconstructed again. Because the 
projection fails with probability 1 — F, the approximate preparation must be 
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done 1/(1 — F) times in average. This can be improved using phase amphfica- 
tion (a trick similar to Grover's search) to ^1/(1 — F) "coherent" initial state 
preparations. A different method is possible if, as is often the case, we can 
evolve with a Hamiltonian H for which the state approximation is a ground 
state. Assume that the approximated ground state has an energy gap bounded 
by A for H and the exact ground state has a similar gap for H. Then we can 
transform a single preparation of the approximated state into the exact ground 
state using around 1/(1 — F) phase estimations, each implemented with a time 
evolution for a time of 1/A (Boixo, Knill & Somma 2010). 

Therefore, a quantum computer, even if it can not solve all ground-state 
problems efficiently, is capable to leverage classical trial states, and solve a 
border class of problems than those efficiently solvable by classical computers. 

Adiabatic state preparation 

The adiabatic method is an alternative way to prepare ground states (Farhi, 
Goldstone, Gutmann & Sipser 2000, Aharonov & Ta-Shma 2003, Perdomo, 
Truncik, Tubert-Brohman, Rose & Aspuru-Guzik 2008, Boixo, Knill & Somma 
2009, Biamonte, Bergholm, Whitfield, Fitzsimons & Aspuru-Guzik 2010, Boixo 
et al. 2010). The original idea is due to Farhi et al. (2000). We first must be able 
to efficiently prepare the ground state |'0 (0)) of a certain initial Hamiltonian 
H{0) — Hi. Then we change the Hamiltonian H{t) slowly, e.g., 

H{t) = {l-t/T)H, + {t/T)Hj. (57) 

Notice that for many reasonable choices of Hi and most physical Hamiltoni- 



2.2 



ans Hf the Hamiltonian H(t) can be simulated using the methods of Sec. 
Nevertheless, common two-body Hamiltonians could be simulated directly. 

If the change from Hi (when i = 0) to the target Hamiltonian Hf (when 
t = T) is slow enough, then the state \ip{t)), satisfying the time-dependent 
Schrodinger equation 

zhj^\i^{t))^H{t)\ij{t)) , (58) 

follows the corresponding eigenstate of H{t) adiabatically. This means that 
\ijj (T)) is close to the ground state of the target Hamiltonian Hf. A sufficient 
condition for the total time T to ensure the adiabaticicty for a linear interpola- 
tion between two Hamiltonians is 

T > ——^ , (59) 

niin 



where s = t/T. Here 



A,„i„= min {E,{s)^Eo{s)) (60) 

0<s<l 



is the minimum gap between the instantaneous eigen-energies Ei{s) and Eo{s) 
of the first excited state and the ground state. The following bound has a 
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better dependence on the minimum gap and it also holds for general (non-linear) 
interpolations if the rate of change of the instantaneous eigenstate \dstp{s)) is 
known (Boixo et al. 2009) 

T>i—- (61) 

Here C is the path length given by the equatiorj^ 

C= f \\dsij{s))\\ds . (62) 



Using the methods of Sec. |2.2| adiabatic evolutions can be simulated effi- 
ciently on a quantum circuit. That is, for cases where one may not be able to 
physically implement H{t), it is still possible to turn the adiabatic state prepa- 
ration into a quantum algorithm and simulate the adiabatic process in a digital 
quantum computer. Furthermore, in this case the total time of the adiabatic 
evolution can be improved tcP] (Boixo et al. 2010): 

T > -^ . (63) 

The remaining question is, in terms of finding ground states, 'how good are 
adiabatic algorithms?'. As we have seen, the performance, or computational 
complexity, of adiabatic algorithms generically depends on the scaling of the 
minimal gap Amjn. Even for classical target Hamiltonians Hf, whether adiabatic 
algorithms success in solving N P-problems is still a controversial issue ( Altshuler, 
Krovi & Roland 2010, Knysh & Smelyanskiy 2010). Numerical results sug- 
gest that for the classical satisfiability (SAT) problems, the scaling of the gap 
would be exponential (Young et al. 2008, Young & Smelyanskiy 2010, Hen & 
Young 2011). If the target Hamiltonian is quantum, the problem is QMA- 
complete. Nevertheless, we can in principle apply the adiabatic algorithm to 
the trial states to improve the ground-state fidelity (Oh 2008), which gives us 
higher probability to project into the exact ground state by the phase estimation 
algorithm discussed in the previous section. 

2.3.2 Preparing thermal states using quantum Metropolis 

We now consider the preparation of density matrices for thermal states 

Pt, = e-^"/TT{e-P") , (64) 

where H can be a quantum or classical Hamiltonian, and /3 = l/T is the inverse 
temperature. We simplify the notation by choosing our units so that the Boltz- 
mann constant fc^ is 1. Zalka (1998b) and Terhal & DiVincenzo (2000) proposed 
to simulate the Markovian dynamics of the system by modeling the interaction 



*More precisely, for this equation we must make a choice of phases such that {dsg{s)\g{s)). 
^Boixo & Somma (2010) have shown that this expression for the total evolution time is 
also optimal. 
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with a heat-bath by some ancilla qubits. A similar idea has been recently in- 
vestigated by Wang, Ashhab & Nori (2011). Terhal & DiVincenzo (2000) also 
attempted to prepare thermal states by generalizing classical Metropolis-type 
sampling (Gould, Tobochnik & Christian 2007). This first quantum Metropolis 
algorithm was limited by the fact that it was not possible to control the up- 
date rule for the Metropolis algorithm, which would generally lead to a slow 
convergence rate of the underlying Markov chain. A significant improvement 
upon this work has been presented recently in Temme, Osborne, Vollbrecht, 
Poulin & Verstraete (2011) with the "quantum Metropolis sampling" algo- 
rithm. This algorithm also constructs a Markov-chain whose fixed point is a 



thermal state Eq. (64 1, but the transitions between states can be engineer to 
achieve faster convergence. The underlying time cost of this algorithm scales 
as 0(1/ A) (Aldous 1982), where A is the eigenvalue gap of the Markov matrix 
associated with the Metropolis algorithm. 

Szegedy (2004) introduced a quantum algorithm to speedup classical Markov 
chains. Richter (2007) extended this methods to some quantum walks with 
decoherence. Szegedy's method has also been applied to achieve a quadratic 
speedup of classical simulated annealing algorithms (Somma, Boixo, Barnum 
& Knill 2008, Wocjan & Abeyesinghe 2008). Further, Yung & Aspuru-Guzik 
(2012) achieved a similar speedup of the quantum Metropolis sampling algo- 
rithm for quantum Hamiltonians. This algorithm outputs a coherent encoding 
of the thermal state (GETS): 



k 

where Z is the partition function, and Ek is the eigenvalue associated with the 
eigenvector jcfc) of the Hamiltonian H. 

Markov chain approaches are practical for many applications. However, for 
systems like spin glasses, the eigenvalue gap A of the corresponding Markov 
matrices typically become exponential small, making it inefficient. Several al- 
ternative approaches have been already introduced in the literature. Among 
them: exploiting the transfer-matrix structure of spin systems (Lidar & Biham 
1997, Yung et al. 2010); mapping the GETS as the ground state of certain 
Hamiltonian (Somma, Batista & Ortiz 2007, Ohzeki & Nishimori 2011); and 
methods based on quantum phase estimation (Poulin & Wocjan 2009b, Bilgin 
& Boixo 2010, Riera, Gogolin & Eisert 2011). In the next section we modified 
one of these algorithms to prepare thermal states building up from small to 
bigger subsystems. 

2.3.3 Preparing thermal states with perturbative updates 

The quantum Metropolis algorithms of the previous subsection extend the ad- 
vantages of Markov chain Monte Garlo methods to the quantum case, even if we 
do not know how to diagonalize the quantum Hamiltonian. It is expected that, 
as in the classical case, they will exhibit good performance for most Hamil- 
tonians. Nevertheless, for very complex systems, such as strongly correlated 
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Hi if 2 H = Hi + if 2 




Hi ^ ii2 H = Hi + H2 



Figure 2: Pictorial representation of the perturbative update method. The 
top figure depicts two quantum systems with Hamihonians iii and ii2 whose 
thermal states we can prepare (maybe through prior perturbative updates). 
The bottom figure depicts a single system where the two quantum systems of 
the top figure have been linked with Hamiltonian h. The perturbative update 
technique is a method to prepare the thermal state of the linked system from 
thermal states of the smaller subsystems. 



molecules, it might be difficult to design rules to choose appropriate Markov 
chain update candidate states, or the convergence rate to the thermal state 
might still be too slow. In this subsection we will show how, even in the worst 
case, quantum algorithms for preparing thermal states will exhibit a substan- 
tial speedup over classical algorithms, elaborating upon the method of Bilgin & 
Boixo (2010). 

The core of this algorithm is a perturbative update subroutine that builds 
the thermal state p'-'^-' ex e"'^'^'*''^''' from the state p'-^' ex e^^^ . We can use 
this subroute to build thermal states of complex systems out of thermal states 
of their components (see Fig. [2|. For this we start with the thermal states 
of subsystems with Hamiltonians Hi and ii2, and use them to prepare the 
thermal state of the linked system with Hamiltonian Hi + H2 + h. The coupling 
Hamiltonian h is introduced perturbatively with a sequence of updates: 

p(o) ^ p(0 ^ p(2.) ^ . . . ^ pH) . (66) 

Quite possibly the thermal states of the smaller subsystems have themselves 
been prepared with a perturbative update over still smaller pieces. As in the 
quantum Metropolis case, it is not necessary to know how to diagonalize the 
corresponding Hamiltonian. 

The perturbative update subroutine is probabilistic, and succeeds with prob- 
ability 1 — e(3Tr p^^^h, which gives the dominant cost of the algorithm. If the 
perturbative update fails, we must reconstruct the state p^*^'. The probability 
of failure is given by the maximum change of a thermal state p*-"' ex e~°^ intro- 
duced by a perturbation eh, which we now bound. We denote with Z = Tr p^*^' 
the partition function of p'-'^'. Using the Dyson series expansion in imaginary 
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time we write 



Z JO 



dX^^-PH(x-x.)^^-pHX,^___ (67) 



The appropriate measure of the difference between two density matrices is the 
trace norm || • Hxr- The reason is that this norm bounds the difference of 
arbitrary measurement results for those two states. The trace norm for an 
arbitrary operator A is given by the sum of the eigenvalues of vMJ\A, and often 
scales with the dimension of the operator. We want to do better for the trace 
norm of the difference between a thermal state and its perturbation, because 
their dimension grows exponentially (in the number of subsystems). We give 
such a bound next. 

We will use the following inequality which applies to all unitarily invariant 
matrix norms l^ 

I A^XB^-^dt <1/2\\\AX + XB\\\ (68) 

Jo 

Applying this inequality to the trace norm of the Dyson series of a perturbed 
thermal state we obtain the bound 

(l/Z) ep f dXie~^"'^^-^''>he-^"^' < eP\\h\\ (69) 

Jo Tr 

where \\h\\ is the operator norm of h. Notice that the operator norm \\h\\ is the 
highest eigenvalue of h, and does not scale with the dimension of H (or even 

The perturbative update subroutine is composed of two operations. The 
first operation implements the quantum map 

(/"'>) ^ (1 - ef3h/2)p^°\l - tPh/2)M , (70) 

where M is just a normalization factor. Similar to the algorithms of the previ- 
ous section, this map is implemented with phase estimation and a conditional 
rotation on an ancillary system. The ancillary system is then measured. This 
measurement can fail, which corresponds to implementing the wrong transfor- 
mation in the thermal state. The success rate is 1 — e/?Trp(°^/i. When the 
measurement of the ancilla system fails, the thermal state can not be recovered, 
and we must start again from the beginning. The cost of the phase estimation 
is e~^/3~^||ft.||~^. This operation can be understood as an update of the Gibbs 
probabilities of p'^^^ to those of p'^^ . The second operation of the perturbative 



I'^See, for instance, Theorem 5.4.7 in (Bhatia 2007). 

^^Although strictly speaking we have derived the bound on the change of the thermal state 
here only to second order, it can be shown to be valid to all orders. For that, we use the exact 
formula for the perturbation e/3 / dX\ e~"^^^~^''-'he~P^^^^i^''- , and the same matrix norm 
inequality. In addition, we need to use the bound for the change on the partition from Poulin 
& Wocjan {2009b). 
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update is a transformation to the eigenbasis of p^'^\ This is implemented by 
"dephasing" in that basis, which is achieved by evolving for a random amount 
of time (with expectation time e||ft.||) using the Hamiltonian H + eh. This com- 
pletes the perturbative update subroutine. 

2.4 Algorithmic quantum cooling 

Yung, Boixo & Aspuru-Guzik (2011) presented an algorithmic quantum cooling 
approach that transforms any input state pi„ into an output state pout which 
is guaranteed to have lower energy with respect to a given Hamiltonian H. 
Explicitly, 

Tr (i?p„„0 < Tr (i/p„) . (71) 

In principle, this algorithm can cool the resulting quantum state to a state 
arbitrarily close to the ground state of H . Nevertheless, like the ground-state 



algorithms of Sec. 2.3.1 the efficiency is related to the energy gap A between 
the ground state and the excited state(s). Depending on how the algorithm is 
implemented, this dependence can scale like 0(1/A^) or 0(1/A). 

Algorithmic quantum cooling first entangles an ancilla qubit with the system 
state. When the ancilla qubit is measured, a result of |0) correlates with a cooler 
system state. On average, however, there is no gain or loss of energy. This 
measurement is used to gain information, just like a Maxwell's demon. The 
measurement outcome of the ancilla qubit in algorithmic quantum cooling can 
be mapped into a ID random walk. The walker starts at a; = 0. For the cooling 
outcome, the walker makes a step towards the positive side x > 0, and towards 
the negative side a; < for the heating outcome. If the walker moves too far 
to the negative side, the procedure is restarted. For some range of parameters, 
whenever the walker goes to the negative side a; < 0, it is guaranteed that the 
quantum state is hotter than the original state. Therefore, removing these hot 
walkers will reduce the average energy over an ensamble of walkers, just like in 
evaporative (or "cofi^ee") cooling of gas molecules. The procedure stops once 
the walker has moved sufficiently to the positive side. 

2.4.1 Basic idea of the quantum cooling method 

We now sketch the basic working mechanism of algorithmic quantum cooling. 
The core component of this cooling algorithm consists of four quantum gates 
(see Fig. [3|p] The first gate is a so-called Hadamard gate: 



H^^(|0) + |1))(0| + ^(|0)-|1))(1| . (72) 

It is followed by a local phase gate, 

R,(7)^|0)(0|-*e'^|l)(l| , (73) 

where the parameter 7 plays a role in determining the overall efficiency of the 
cooling performance of the algorithm. The interaction with the Hamiltonian 



^^Similar quantum circuits arc used in DQCl and phase estimation, for instance. 
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Ancilla |0)- 



Input state 



Rz(v) 



Measurement 



Output state 



Figure 3: Quantum circuit diagram for the quantum cooling algorithm. Here 



H 



)= (|0) + |1)) (0| + ^ (|0) - |1)) (1| is the Hadamard gate, Rz (7) = |0) 



V2 
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|1)(1| is a local phase gate, and U {t) — 



-,-iH^t 



, is the time t evolution 



operator simulating the dynamics of the system under the Hamiltonian Hg ■ 

H, which can be either quantum or classical, is encoded in the time evolution 
operator 

U (t) = e-*^=* . (74) 

As explained above, time evolution can be implemented efficiently in a quantum 
computer. 

The operation of the circuit in Fig. Islon input state |V'm) is as follows. 

Step 1 State initialization, 

\An) |0) , (75) 

with the ancilla state in |0). 

Step 2 Apply the Hadamard gate, H = ^ (|0) + |1)) (0| + ^ (|0) - |1)) (1|, 
and the local phase gate R^ (7) = |0) (0| — ie^'' |1) (1| to the ancilla qubit, 

|Vm)(|0)-*e*^|l))/%/2 . (76) 

Step 3 Apply the controlled-t/(t) to the system state, 

(|i^„,)|0)-ze^'^[/(i)|V™)|l))/\/2 . (77) 



Step 4 Apply the Hadamard to the ancilla qubit again. This produces the 
following output state: 



Ao|V»n)|0)+Ai|V^™)|l) , 



(78) 



where Aj = (/ + {-iy+^ie''^U) /2 for j = {0, 1}. 

A projective measurement on the ancilla qubit in the computational basis 
{|0) , |1)} yields one of the two (unnormalized) states. 



(/±ze^-^C/)|V„) 



(79) 



Their mean energy is either higher (for outcome |1), x is decreased by 1) or 
lower (for outcome |0), x is increased by 1) than that of the initial state \ipm)- 
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To justify this assertion, let us expand the input state, 

\ipin) = V", Cfe \ek) , (80) 

in the eigenvector basis {|efe)} of the Hamiltonian H. Note that 

|(l±ze*TC/)|efe)|' = |cfe|'(l±sin(/)fc) , (81) 

where 4'k = Ekt — 7 depends on the eigen-energy Ek of H. For simplicity, we 
will assume that one can always adjust the two parameters, 7 and t, such that 

- f < .^fc < f (82) 

for all non-negative integers k. Then, the factors (1 — siiKpk) are in descending 
order of the eigen-energies, and the opposite is true for the factors (1 + sint/)^). 
Therefore, apart from an overall normalization constant, the action of the opera- 
tor (/ ± ie^^U) is to scale each of the probability weights \ck\ by an eigen-energy 
dependent factor (1 ± sin^^), i.e., 

|cfe|'^|cfc|'(l±sin,/)fe) . (83) 

The probability weights scale to larger values, i.e., 

(l-sin(/)fc)/(l-sin0j) > 1 (84) 

for the eigen-energy Ek < Ej in the cooling case (i.e., for outcome |0)), and 
vice versa for the heating case (i.e., for outcome |1)). Further cooling can 
be achieved by applying the quantum circuit repeatedly and reject/recycle the 
random walker when x < 0. 

2.4.2 Connection with heat-bath algorithmic cooling 

The algorithmic quantum cooling approach is related to the well-known heat- 
bath algorithmic cooling (HBAC) (Boykin, Mor, Roychowdhury, Vatan & Vrijen 
2002, Baugh, Moussa, Ryan, Nayak & Laflamme 2005, Schulman, Mor & Weinstein 
2005). HBAC aims to polarize groups of spins as much as possible, i.e. to pre- 
pare the state 

\m - 1) . (85) 

This state is important for providing fresh ancilla qubits for quantum error 
correction as well as for NMR quantum computation. In HBAC, some reversible 
operations are first performed to redistribute the entropy among a group of 
spins. Some of the spins will become more polarized. For a closed system, 
there is a so-called Shannon bound (Schulman et al. 2005) which limits the 
compression of the entropy. In order to decrease the entropy of the whole system, 
the depolarized spins interact with a physical heat bath that acts as an entropy 
sink. We note that from an algorithm point of view, the existence of a physical 
heat bath can be replaced by the (imperfect) preparation of polarized spins by 
other methods. The method of algorithmic quantum cooling from Yung et al. 
(2011) may be considered as a generalization of the HBAC, as it is applicable 
to cool any physical system that is simulable by a quantum computer, not just 
non-interacting spins. 
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3 Special topics 

3.1 Adiabatic non-destructive measurements 



In Section. 2.3 we reviewed several methods to prepare ground states and ther- 
mal states of quantum systems of interest in physics and chemistry. In par- 
ticular, in subsection |2.3.1| we gave an overview of the adiabatic method for 
preparing ground states. The adiabatic model may be naturally more robust 
against noise, offering a method to perform small to medium size simulations 
without using sophisticated error correction schemes. Because of this and other 
reasons, adiabatic based quantum computation is possibly easier to realize phys- 
ically than quantum computation based on the circuit model. In this section 
we review a method to effect non-destructive measurements of constants of the 
motion within the adiabatic model. 

As explained in subsection |2.3.1[ it is in principle possible to adiabatically 
prepare the ground state of a physical or chemical system with Hamiltonian 
Hf. There we said that this can be done by interpolating slowly enough be- 
tween a simple initial Hamiltonian Hi and the final Hamiltonian Hf . Following 
Biamonte et al. (2010), we now add an ancillary qubit subsystem with orthonor- 
mal basis {\po), \pi)}- This auxiliary system will be use for the adiabatic non- 
destructive measurements. During the adiabatic ground state preparation, this 
subsystem is acted upon with Hamiltonian i5|pi)(pi|, and therefore it remains 
in the state \po)- The choice of i5 > will be explained shortly. 

The measurement procedure begins by bringing the ancillary qubit and the 
system being simulated into interaction, adiabaticalljj^ We choose the interac- 
tion Hamiltonian Hint = ^'^ \Pi){pi\- Here A is any observable corresponding 
to a constant of the motion, that is [A,H] = 0. In particular, the Hamiltonian 
Hf itself can be used to obtain the ground state energy. The total Hamiltonian 
becomes 

Hf + 6\pi){pi\+A®\pi){pi\ . (86) 

Hsp 

If the energy bias S is bigger than the expectation value of the observable A, 
the state does not change during this initial interaction (Biamonte et al. 2010). 
After the initial interaction, we apply a Hadamard gate to the ancillary 
qubit. We denote the time at which we apply this gate as i = 0. Let |so) be the 
ground state of Hf. After a further time t the system plus ancilla qubit evolves 
to 

|^(t)) = ^|so)®(bo)+e-'-*bi)) (87) 

where uj = (a^ + S)/h, and oq = (sqI^Iso) is the expectation value we wish to 
measure. Finally, we again apply a Hadamard gate to the probe. The resulting 



^^The interaction Hamiltonian is typically a three-body Hamiltonian, which makes direct 
simulations more difficult. This difficulty can be overcome using gadgets (Kempe, Kitaev & 
Regev 2006b, Oliveira & Terhal 2005, Biamonte & Love 2008, Biamonte et al. 2010) or the 
average Hamiltonian method (Waugh, Huber & Haeberlen 1968) 
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state is 

\^it)) = |so> ® (cos {ujt/2) Ipo) + i sin (cjt/2) |pi)) , (88) 

yielding probability, 

Po{t) == o (1 + cos(a;i)) = cos^ (cjt/2) . (89) 

Measuring the probe does not disturb the state of the simulator which can 
be reused for another measurement. This measurement can be repeated until 
sufficient statistics have been accumulated to reconstruct to. We refer to Bia- 
monte et al. (2010) for details on numerical simulations and considerations of 
the influence of noise. 

3.2 TDDFT and quantum simulation 

Density Functional Theory (DFT) and its time-dependent extension (TDDFT) 
have become arguably the most widely used methods in computational chem- 
istry and physics. In DFT and TDDFT, the properties of a many-body system 
can be obtained as functionals of the simple one-electron density rather than the 
correlated many-electron wavefunction. This represents a great conceptual leap 
from usual wavefunction-based methods such as Hartree-Fock, configuration in- 
teraction and coupled cluster methods and therefore the connections between 
DFT/TDDFT and quantum computation have just begun to be explored. Since 
TDDFT is a time-dependent theory, it is more readily applicable to quantum 
simulation than DFT which is strictly a ground state theory. For recent devel- 
opments in the connections between DFT and quantum complexitjj^see Schuch 
& Verstraete (2009b), while for applications of DFT to adiabatic quantum com- 
putation see Gaitan & Nori (2009). In this section, we provide a brief overview 
of the fundamental theorems of TDDFT, which establish its use as a tool for 
simulating quantum many-electron atomic, molecular and solid-state systems 
and we mention recent extensions of TDDFT to quantum computation (Tempel 
& Aspuru-Guzik 2011). 

In its usual formulation, TDDFT is applied to a system of N-electrons de- 
scribed by the Hamiltonian 

W -2 ^ /. 

^W = E 2^ + E^(|f^ - ^^1) + / v{r,t)niv)dh, (90) 

where Pi and f^ are respectively the position and momentum operators of 
the jth electron, w(|fi — fj|) is the electron-electron repulsion and w(r,t) is 
a time-dependent one-body scalar potential which includes the potential due 
to nuclear charges as well as any external fields. The electron-electron repul- 
sion, w(|fi — fjl), leads to an exponential scaling of the Hilbert space with 
system-size and makes simulation of the many-electron Schrodinger equation 



^*It turns out that finding a universal funcional for DFT is QMA hard. Liu, Christandl 
Verstraete (2007a) proved a related results: A'-representability is also QMA hard. 
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on a classical computer intractable. n(r) = J2i H'"^ ~ ^i) is the electron density 
operator, whose expectation value yields the one-electron probability density, 
(n(r)) = n{r,t). The basic theorems of TDDFT prove that, in principle, one 



can simulate the evolution of the Hamiltonian in Eq. 90 using n(r, t) directly and 
thereby avoid calculating and storing the exponential amount of information in 
the many-electron wavefunction. 

The first basic theorem of TDDFT, known as the "Runge-Gross" (RG) the- 
orem" (Runge & Gross 1984), establishes the existence of a one-to-one mapping 
between the expectation value of n(r) and the scalar potential v{r,t). i.e 

n{r,t)^v{r,t). (91) 



However, v{r,t) is the only part of the Hamiltonian in Eq. 90 that is non- 



universal, i.e. J2i=i 2m + Si<i ^(l^j ~ ^j\) is the same operator for each 
electronic system. Therefore, due to the uniqueness of the solution to the 
time-dependent Schrodinger equation, the RG theorem establishes a one-to- 
one mapping between the density and the wavefunction. This implies that the 
wavefunction is in fact a unique functional of the density, 

V'(ri,---,rAr;i) = V'W(ri,...,rAr;t), (92) 

as is any observable of the system. The RG theorem implies the remarkable 
fact that the one-electron density contains the same quantum information as 
the many-electron wavefunction. This means that in principle, if one had a 
means of directly simulating n{r,t), one could extract all observables of the 
system without ever needing to simulate the many-body wavefunction. 

The second basic TDDFT theorem is known as the "van Leeuwen (VL) the- 
orem" (van Leeuwen 1999). It gives an analytic expression for a time-dependent 
one-body scalar potential that applied to another system with a different, and 
possibly simpler, electron-electron repulsion ui'df, — fj|), gives the same density 
evolution as the original Hamiltonian of Eq. 90 When w'{\ri—rj\) = 0, this aux- 
iliary system is referred to as the "Kohn-Sham system" (Kohn & Sham 1965). 
Due to it's simplicity and accuracy, the Kohn-Sham system is in practice used 
in most DFT and TDDFT calculations. Since the Kohn-Sham system is non- 
interacting, its wavefunction is simply described by a Slater determinant of 
single-electron orbitals, which satisfy the time-dependent Kohn-Sham equations. 






-7;'^^+Vks[n]{r,t) 



Mr,t). (93) 



The true interacting density is obtained from the orbitals by square-summing; 
that is, n{r,t) = J2i l'/'i(i"jOP- Naturally, the set of single-particle equations in 
Eq. [93] are far easier to solve than evolution under the Hamiltonian in Eq. [90) 
In practice, the Kohn-Sham potential, Vks[n]{r,t), must be approximated as a 
density functional, but the VL theorem rigorously guarantees its existence and 
uniqueness. 
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1-1 mapping: RG Theorem 



b) lH[ai,ala^{t) 



/i2[o-J,(T|,cr|](t) 



h3lo't,o-l,o-i]{t) 






I 



1-1 mapping: Schrodinger Equation 



|^(t)) = |V^[af,a|,af](t)) 



Figure 4: Runge-Gross theorem for a 3 qubit example - The set of ex- 
pectation values {erf ,ct|, ...ct^}, defined by the the Bloch vector components of 
each qubit along the z-axis in (a), is uniquely mapped onto the set of local fields 
{/ii,/i2, ...ft-Af} in (b) through the RG theorem. Then, through the Schrodinger 
equation, the set of fields is uniquely mapped onto the wavefunction. These 
two mappings together imply that the N-qubit wavefunction in (c) is in fact a 
unique functional of the set of expectation values {erf , cr|, ...cr^}. 

Tempel & Aspuru-Guzik (2011) recently extended the RG and VL theo- 
rems to systems of interacting qubits described by the class of universal 2-local 
Hamiltonians 



h,it)a^. (94) 



This Hamiltonians apply to a variety of different systems, particularly in solid- 
state quantum computing. Benjamin & Bose (2003) and Benjamin & Bose 
(2004) have shown that any set of universal two-qubit and single-qubit gates 
can be implemented with the Hamiltonian of Eq. |94[ and therefore it can be 
used to perform universal quantum computation. The RG theorem applied to 
such universal Hamiltonians establishes a one-to-one mapping between the set 





JV-l W-l N 


H{t). 


- E '^+i(-r-?+i + ^r-^i) + E ^!.+i^r-^i + E 
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J 1-1 mapping: VL Theorem {^^ ^4'} + {Ky. 4^ 



b) 



''i[o-r,o-5,<T5](t) 




UJya\,al,al\{t) 



h'jyal.al.al\[t) 



cz^ c±:> 



J 1-1 mapping: Schrodinger Equation |-i^'(i)^ ^ | j_;i(f)^ 




^2~ = f^2 



^Z = (^1 





Figure 5: Van Leeuwen theorem for a 3 qubit 

{erf ,(t|, ■■■<jff} (a) obtained from evolution under Eq.|94 



example - The set 

is uniquely mapped to 
a new set of fields {h'l, h'2, ■■■h'j^} (b) for a Hamiltonian with different two-qubit 
interactions. Evolution under this new Hamiltonian returns the same expec- 
tation values {erf , erf, ...CT^}, although the wavefunction is different and hence 
projections of the Bloch vectors along other axes are in general different (c). 



of local fields {hi, h2, .../iat} used to implement a given computation, and the set 
of single-qubit expectation values {cf ,(t|, ••■<t^} (see Fig. 111). This implies that 
one can use single-qubit expectation values as the basic variables in quantum 
computations rather than wavefunctions and directly extract all observables of 
interest with only knowledge of the spin densities. Naturally, certain properties 
such as entanglement will be difficult to extract as functionals of the set of spin 
densities {crf,cr|, ...(jf^}. Nevertheless, Tempel & Aspuru-Guzik (2011) give an 
explicit entanglement functional. 

In addition to the RG theorem, one can derive a VL theorem for qubits. 
The VL theorem for qubits provides an exact prescription for simulating uni- 
versal Hamiltonians with other universal Hamiltonians that have different, and 
possibly easier-to-realize, two-qubit interactions. In analogy to the Kohn-Sham 
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system in electronic TDDFT, one can consider an auxiliary Hamiltonian, 

H'(t) = E'^^i('^?'^^i + '^'^m) + E'^l+i'^^^m 

i=l i=l 

N 

+ Y.^[K,al,...a^^]{t)al (95) 

i=l 

with simpler two-qubit couplings {J'-*", J'"}. The VL theorem guarantees the 
existence of the auxiliary fields {h'i,h'2, ■■■h'j^} as functional of the spin densities 
which reproduce any given set {af, (t|, •■■<t^} that one might wish to simulate. 
In this way, one can construct an entire class of density functionals for quantum 
computing that map between different universal Hamiltonians as illustrated in 
figure [5] 

TDDFT applied to universal qubit Hamiltonians provides a potentially pow- 
erful tool for simulating quantum computations on classical computers, similar 
to how it has been applied in computational chemistry for simulating electronic 
systems. By choosing the auxiliary "Kohn-Sham" system to be less entan- 
gled than the original system, one can hope to simplify simulations of quan- 
tum algorithms by finding simple approximations to the auxiliary local fields 
{h[, /12, ■■■h'j^} as functionals of the spin density. As in electronic TDDFT, the 
development of approximate density functionals for qubit systems will be a nec- 
essary next step, which is discussed in Tempel & Aspuru-Guzik (2011). 

4 Conclusion and outlook 

To the best of our knowledge, the first quantum simulation experiment was 
performed by Somaroo, Tseng, Havel, Laflamme & Cory (1999) in a two-qubit 
NMR system in 1999, where a truncated quantum harmonic or anharmonic os- 
cillator (4 levels) was simulated. Strictly speaking, this belongs to an "analog" 
simulation; because the Hamiltonian of the quantum oscillator was directly sim- 
ulated by the Hamiltonian of the spins, instead of applying quantum algorithms. 
The progress of quantum simulation is still gaining momentum. In Table [ll we 
list several recent experiments on digital quantum simulation. Earlier references 
may be found from them and Buluta & Nori (2009), Kassal et al. (2011), and 
Jones (2011). 

So far, none of the simulations implement active error correction. An im- 
portant aspect of digital quantum simulation is the resource estimation when 
fault-tolerant structures are considered. We refer to Brown, Clark & Chuang 
(2006) and Clark et al. (2009) for further explore this area. In short, to achieve 
large-scale quantum simulation, there are still many technological challenges to 
overcome. For example, in many currently available setups, the performances 
of the two-qubit gates are still too noisy for fault-tolerant simulation. However, 
by experimenting with small-scale quantum simulation experiments, we believe 
that valuable lessons can be learnt to optimize the performance of the currently 
available technology. 
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Table 1: A brief survey of recent experiments on digital quantum simulation 

Physical implementations 



What is simulated? 



Scale 



Nuclear Magnetic Resonance (NMR) 



Trapped Ions 



• Thermal states of a frustrated magnet (Zhang, 4 qubits 
Yung, Laflamme, Aspuru-Guzik h Baugh 2011) 

• Ground state of a pair of interacting Hciscn- 3 qubits 
berg spins subject to simulated external fields (Li, 

Yung, Chen, Lu, Whitfield, Peng, Aspuru-Guzik 
k Du 2011) 

• Isomcrization reaction dynamics (Lu, Xu, Xu, 3 qubits 
Chen, Gong, Peng & Du 2011) 

• Ground state of hydrogen molecule H2 (Du, Xu, 2 qubits 
Peng, Wang, Wu & Lu 2010) 

• Time dynamics of spin systems (Lanyon, 6 qubits 
Hempel, Nigg, MuUer, Gerritsma, Zahringer, 
Schindler, Barreiro, Rambach, Kirchmair, Hcn- 

nrich, ZoUer, Blatt k Roos 2011) 

• Dissipative open-system dynamics (Barreiro, 5 qubits 
Miiller, Schindler, Nigg, Monz, Chwalla, Hcnnrich, 

Roos, ZoUer k Blatt 2011) 



Quantum Optics 



• ID quantum walk of a topological sys- 4 steps 
tem (Kitagawa, Broome, Fedrizzi, Rudner, Berg, 
Kassal, Aspuru-Guzik, Demlcr k White 2011) 

• ID quantum walk with tunable decoher- 6 steps 
ence (Broome, Fedrizzi, Lanyon, Kassal, Aspuru- 
Guzik k White 2010) 

• Ground state of hydrogen molecule H2 (Lanyon, 2 qubits 
Whitfield, Gillett, Goggin, Almeida, Kassal, Bia- 

montc, Mohscni, Powell, Barbieri, Aspuru-Guzik 
k White 2010) 



A related question is "what is the best way to implement quantum simu- 
lation?". For classical computers, there is no doubt that silicon-based semi- 
conductors work successfully. For quantum computers, a general feature of the 
currently proposed technologies, such as quantum dots, quantum optics, trapped 
ions, nuclear and electron spins, impurity, superconducting devices etc., is that 
there is a trade-off between controllability and reliability. Usually, systems that 
can be controlled easily suffer more from decohcrence from the environment. 
There are two approaches to tackle this problem: one may cither look for new 
systems that are good for both control and can be isolated from the environment, 
or develop hybrid structures that combine the advantages from both sides. For 
example, there has been progress in coupling superconducting devices with spin 
ensembles (Duty 2010). The former provides the controllability and the latter 
provides reliability. In short, the future of quantum computation and quantum 
simulation is still full of challenges and opportunities. We hope this article can 
stimulate more ideas that can help move this field forward. 



34 



5 Acknowledgements 

The authors would hke to acknowledge NSF CCI grant number: CHE-1037992. 

References 

Aaronson, S. (2009). Computational complexity: Why quantum chemistry is 
hard, Nature Physics 5(10): 707-708. 

Aaronson, S. (2011). Why Philosophers Should Care About Computational 
Complexity, arXiv:1108.1791 p. 58. 

Abrams, D. & Lloyd, S. (1997). Simulation of many-body Fermi systems on a 
universal quantum computer, Phys. Rev. Lett. 79: 2586. 

Abrams, D. & Lloyd, S. (1999). Quantum Algorithm Providing Exponential 
Speed Increase for Finding Eigenvalues and Eigenvectors, Physical Review 
Letters 83(24): 5162-5165. 

Aharonov, D. & Naveh, T. (2002). Quantum NP - A Survey, arXiv.-quant- 
ph/0210077p. 23. 

Aharonov, D. & Ta-Shma, A. (2003). Adiabatic quantum state generation and 
statistical zero knowledge. Proceedings of the thirty-fifth annual ACM sym- 
posium on Theory of computing, ACM, pp. 20-29. 

Aldous, D. J. (1982). Some Inequalities for Reversible Markov Chains, Journal 
of the London Mathematical Society s2-25(3): 564-576. 

Altshuler, B., Krovi, H. & Roland, J. (2010). Anderson localization makes adia- 
batic quantum optimization fail, Proc Natl Acad Sci USA 107(28): 12446- 
50. 

Aspuru-Guzik, A., Dutoi, A. D., Love, P. J. & Head-Gordon, M. (2005). Sim- 
ulated quantum computation of molecular energies.. Science (New York, 
N.Y.) 309(5741): 1704-7. 

Bacon, D., Childs, A., Chuang, I., Kempe, J., Leung, D. & Zhou, X. (2001). 
Universal simulation of Markovian quantum dynamics. Physical Review A 
64(6). 

Barahona, F. (1982). On the computational complexity of ising spin glass mod- 
els, J. Phys. A: Math. Gen 15(10): 3241-3253. 

Barreiro, J. T., Miiller, M., Schindler, P., Nigg, D., Monz, T., Chwalla, M., 
Hennrich, M., Roos, C. F., ZoUer, P. & Blatt, R. (2011). An open-system 
quantum simulator with trapped ions, Nature 470(7335): 486-491. 



35 



Baugh, J., Moussa, O., Ryan, C. a., Nayak, A. & Laflamme, R. (2005). Exper- 
imental implementation of heat-bath algorithmic cooling using solid-state 
nuclear magnetic resonance., Nature 438(7067): 470-3. 

Deals, R., Buhrman, H., Cleve, R., Mosca, M. & de Wolf, R. (1998). Quantum 
lower bounds by polynomials. Proceedings of FOCS' 98, pp. 352-361. 

Benenti, G. & Strini, G. (2008). Quantum simulation of the single-particle 
Schrodinger equation, American Journal of Physics 76(7): 657. 

Benjamin, S. G. & Bose, S. (2003). Quantum computing with an always-on 
heisenberg interaction, Phys. Rev. Lett. 90: 247901. 

Benjamin, S. G. & Bose, S. (2004). Quantum computing in arrays coupled by 
"always-on" interactions, Phys. Rev. A 70: 032314. 

Berry, D. W., Ahokas, G., Gleve, R. & Sanders, B. G. (2006). Efficient Quan- 
tum Algorithms for Simulating Sparse Hamiltonians, Communications in 
Mathematical Physics 270(2): 359-371. 

Berry, D. W. & Ghilds, A. M. (2009). Black-box Hamiltonian simulations and 
unitary implementation, arXiv:0910.4157 . 

Bhatia, R. (2007). Positive definite matrices, Princeton University Press. 

Biamontc, J. D., Bergholm, V., Whitfield, J. D., Fitzsimons, J. & Aspuru-Guzik, 
A. (2010). Adiabatic Quantum Simulators, University Computing p. 9. 

Biamonte, J. D. & Love, P. J. (2008). Realizable hamiltonians for universal 
adiabatic quantum computers, Physical Review A (Atomic, Molecular, and 
Optical Physics) 78(1): 012352-7. 

Bilgin, E. & Boixo, S. (2010). Preparing Thermal States of Quantum Systems 
by Dimension Reduction, Physical Review Letters 105(17). 

Boixo, S., Knill, E. & Somma, R. D. (2009). Eigenpath traversal by phase 
randomization. Quantum Information and Computation 9: 833-855. 

Boixo, S., Knill, E. & Somma, R. D. (2010). Fast quantum algorithms for 
traversing paths of eigenstates, arXiv:1005.3034 ■ 

Boixo, S. & Somma, R. D. (2010). Necessary condition for the quantum adia- 
batic approximation. Physical Review A 81(3): 032308. 

Boykin, P. O., Mor, T., Roychowdhury, V., Vatan, F. & Vrijen, R. (2002). Algo- 
rithmic cooling and scalable NMR quantum computers.. Proceedings of the 
National Academy of Sciences of the United States of America 99(6): 3388- 
93. 

Breuer, H.-P. & Petruccione, F. (2002). The theory of open quantum systems, 
Oxford University Press. 



36 



Broome, M. A., Fcdrizzi, A., Lanyon, B. P., Kassal, I., Aspuru-Guzik, A. & 
White, A. G. (2010). Discrete Single-Photon Quantum Walks with Tunable 
Decoherence, Physical Review Letters 104(15): 1-4. 

Brown, K. L., Munro, W. J. & Kcndon, V. M. (2010). Using Quantum Com- 
puters for Quantum Simulation, Entropy 12(11): 2268-2307. 

Brown, K. R., Clark, R. J. & Chuang, I. L. (2006). Limitations of quantum 
simulation examined by simulating a pairing Hamiltonian using nuclear 
magnetic resonance.. Physical Review Letters 97(5): 050504. 

Buluta, I. & Nori, F. (2009). Quantum simulators.. Science (New York, N.Y.) 
326(5949): 108- 11. 

Byrnes, T. & Yamamoto, Y. (2006). Simulating lattice gauge theories on a 
quantum computer. Physical Review A 73(2): 1-16. 

Childs, A. M. (2004). Quantum information processing in continuous time, 
Ph.D., MIT, Cambridge, MA. 

Childs, A. M. (2010). On the relationship between continuous- and discrete-time 
quantum walk, Comm. Math. Phys. 294: 581. 

Childs, A. M. & Kothari, R. (2011). Simulating sparse Hamiltonians with 
star decompositions. Theory of Quantum Computation Communication 
and Cryptography TQC 2010 6519: 94-103. 

Cho, J., Bosc, S. & Kim, M. (2011). Optical Pumping into Many-Body Entan- 
glement, Physical Review Letters 106(2): 1~4. 

Clark, C, Metodi, T., Gasster, S. & Brown, K. (2009). Resource requirements 
for fault-tolerant quantum simulation: The ground state of the transverse 
Ising model. Physical Review A 79(6): 1-9. 

Du, J., Xu, N., Peng, X., Wang, P., Wu, S. & Lu, D. (2010). NMR Implemen- 
tation of a Molecular Hydrogen Quantum Simulation with Adiabatic State 
Preparation, Physical Review Letters 104(3): 1-4. 

Duty, T. (2010). Towards superconductor-spin ensemble hybrid quantum sys- 
tems. Physics 3: 80. 

Farhi, E., Goldstone, J., Gutmann, S. & Sipser, M. (2000). Quantum Compu- 
tation by Adiabatic Evolution, arXiv:quant-ph/0001106vl . 

Feit, M., Fleck Jr, J. & Stciger, A. (1982). Solution of the Schrodinger equation 
by a spectral method. Journal of Computational Physics 47(3): 412-433. 

Feynman, R. P. (1982). Simulating physics with computers. International Jour- 
nal of Theoretical Physics 21(6-7): 467-488. 



37 



Gaitan, F. & Nori, F. (2009). Density functional theory and quantum compu- 
tation, Phys. Rev. B 79: 205117. 

Gorini, V., Kossakokowski, A. & Sudarshan, E. C. G. (1976). Completely posi- 
tive dynamical semigroups of n-level systems, J. Math. Phys 17: 821. 

Gould, H., Tobochnik, J. & Christian, W. (2007). An introduction to com- 
puter simulation methods: applications to physical systems, Pearson Addi- 
son Wesley. 

Grover, L. & Rudolph, T. (2002). Creating superpositions that correspond to 
efficiently integrable probability distributions, arXiv:quant-ph/0208 11 2 . 

Hatano, N. & Suzuki, M. (2005). Finding exponential product formulas of higher 
orders, in A. Das & B. Chakrabarti (eds). Quantum Annealing and Other 
Optimization Methods, Lectures Notes in Physics, Springer, pp. 37-68. 

Hauke, P., Cucchictti, F. M., Tagliacozzo, L., Deutsch, I. & Lewenstein, M. 
(2011). Can One Trust Quantum Simulators?, arXiv: 1 109.6457 . 

Helgaker, T., Jorgensen, P. & Olsen, J. (2000). Molecular Electronic- Structure 
Theory, John Wiley and Sons. 

Hen, I. & Young, A. P. (2011). Exponential complexity of the quantum adi- 
abatic algorithm for certain satisfiability problems. Physical Review E 
84(6): 061152. 

Jaksch, P. & Papageorgiou, A. (2003). Eigenvector Approximation Leading to 
Exponential Speedup of Quantum Eigenvalue Calculation, Physical Review 
Letters 91(25): 1-4. 

Jones, J. A. (2011). Quantum computing with NMR., Progress in nuclear mag- 
netic resonance spectroscopy 59(2): 91-120. 

Kassal, L & Aspuru-Guzik, A. (2009). Quantum algorithm for molecular 
properties and geometry optimization.. The Journal of chemical physics 
131(22): 224102. 

Kassal, L, Jordan, S. P., Love, P. J., Mohseni, M. & Aspuru-Guzik, A. (2008). 
Polynomial-time quantum algorithm for the simulation of chemical dynam- 
ics.. Proceedings of the National Academy of Sciences of the United States 
o/ America 105(48): 18681-6. 

Kassal, I., Whitfield, J. D., Perdomo-Ortiz, A., Yung, M.-H. & Aspuru-Guzik, 
A. (2011). Simulating chemistry using quantum computers.. Annual review 
of physical chemistry 62: 185-207. 

Kaye, P., Laflamme, R. & Mosca, M. (2007). An Introduction to Quantum 
Computing, 1 edn, Oxford University Press, USA. 



38 



Kaye, P. & Mosca, M. (2004). Quantum networks for generating arbitrary 
quantum states. 

Kempe, J., Kitaev, A. & Regev, O. (2006a). The complexity of the local hamil- 
tonian problem, SI AM Journal on Computing 35(5): 1070-1097. 

Kempe, J., Kitaev, A. & Regev, O. (2006b). The complexity of the local hamil- 
tonian problem, SI AM Journal on Computing 35(5): 1070. 

Kitaev, A., Shen, A. & Vyalyi, M. (2002). Classical and quantum computation, 
Graduate studies in mathematics, American Mathematical Society. 

Kitagawa, T., Broome, M. A., Fedrizzi, A., Rudner, M. S., Berg, E., Kassal, 
I., Aspuru-Guzik, A., Demler, E. & White, A. G. (2011). Observation of 
topologically protected bound states in a one dimensional photonic system, 
arXiv:1105.5334 p. 5. 

Kliesch, M., Barthel, T., Gogolin, C., Kastoryano, M. & Eisert, J. (2011). 
Dissipative Quantum Church- Turing Theorem, Physical Review Letters 
107(12): 1-5. 

Knill, E., Ortiz, G. & Somma, R. D. (2007). Optimal quantum measurements of 
expectation values of observables. Physical Review A (Atomic, Molecular, 
and Optical Physics) 75(1): 012328-13. 

Knysh, S. & Smelyanskiy, V. (2010). On the relevance of avoided crossings 
away from quantum critical point to the complexity of quantum adiabatic 
algorithm, arXiv:1005.3011 p. 8. 

Kohn, W. (1999). Nobel Lecture: Electronic structure of matter — wave func- 
tions and density functionals. Reviews of Modern Physics 71(5): 1253-1266. 

Kohn, W. & Sham, L. J. (1965). Self-consistent equations including exchange 
and correlation effects, Phys. Rev. 140: 1133. 

Kosloff, R. (1988). Time-dependent quantum-mechanical methods for molecular 
dynamics. The Journal of Physical Chemistry 92{8): 2087-2100. 

Kraus, B., Biichler, H., Diehl, S., Kantian, A., Micheh, A. & ZoUer, P. (2008). 
Preparation of entangled states by quantum Markov processes. Physical 
Review A 78(4): 1-9. 

Krauter, H., Muschik, C., Jensen, K., Wasilewski, W., Petersen, J., Cirac, J. 
& Polzik, E. (2011). Entanglement Generated by Dissipation and Steady 
State Entanglement of Two Macroscopic Objects, Physical Review Letters 
107(8): 1-5. 

Ladd, T. D., Jelezko, F., Laflamme, R., Nakamura, Y., Monroe, C. & O'Brien, 
J. L. (2010). Quantum computers.. Nature 464(7285): 45-53. 



39 



Lanyon, B. P., Hempel, C, Nigg, D., MuUcr, M., Gerritsma, R., Zahringer, 
F., Schindler, P., Barreiro, J. T., Rambach, M., Kirchmair, G., Hennrich, 
M., Zollcr, P., Blatt, R. & Roos, C. F. (2011). Universal Digital Quantum 
Simulation with Trapped Ions, Science p. 13. 

Lanyon, B. P., Whitfield, J. D., Gillett, G. G., Goggin, M. E., Almeida, M. P., 
Kassal, I., Biamonte, J. D., Mohseni, M., Powell, B. J., Barbieri, M., 
Aspuru-Guzik, A. & White, A. G. (2010). Towards quantum chemistry 
on a quantum computer. Nature Chemistry 2(2): 106-111. 

Li, Z., Yung, M.-H., Chen, H., Lu, D., Whitfield, J. D., Peng, X., Aspuru- 
Guzik, A. & Du, J. (2011). Solving Quantum Ground-State Problems with 
Nuclear Magnetic Resonance, Scientific Reports 1: 88. 

Lidar, D. & Biham, O. (1997). Simulating Ising spin glasses on a quantum 
computer. Physical Review E 56(3): 3661-3681. 

Lidar, D. & Wang, H. (1999). Calculating the thermal rate constant with expo- 
nential speedup on a quantum computer. Physical Review E 59(2): 2429- 
2438. 

Lindblad, G. (1975). On the generators of quantum dynamical semigroups, 
Commun. Math. Phys 48: 119. 

Liu, Y., Christandl, M. & Verstraete, F. (2007a). Quantum computational 
complexity of the N-Representability problem: QMA complete, Physical 
Review Letters 98(11): 110503. 

Liu, Y.-K., Christandl, M. & Verstraete, F. (2007b). Quantum Computational 
Complexity of the N-Representability Problem: QMA Complete, Physical 
Review Letters 98(11): 110503. 

Lloyd, S. (1996). Universal Quantum Simulators, S'cience 273(5278): 1073-1078. 

Lu, D., Xu, N., Xu, R., Chen, H., Gong, J., Peng, X. & Du, J. (2011). Simu- 
lation of Chemical Isomerization Reaction Dynamics on a NMR Quantum 
Simulator, Physical Review Letters 107(2): 8-11. 

Master, C, Yamaguchi, F. & Yamamoto, Y. (2003). Efficiency of free-energy 
calculations of spin lattices by spectral quantum algorithms. Physical Re- 
view A 67(3). 

Mohseni, M., Rebentrost, P., Lloyd, S. & Aspuru-Guzik, A. (2008). 
Environment-assisted quantum walks in energy transfer of photosynthetic 
complexes. Journal of Chemical Physics 129(174106). 

Miiller, M., Hammerer, K., Zhou, Y. L., Roos, C. F. & Zoller, P. (2011). Sim- 
ulating open quantum systems: from many-body interactions to stabilizer 
pumping. New Journal of Physics 13(8): 085007. 



40 



Muschik, C, Polzik, E. & Cirac, J. (2011). Dissipatively driven entanglement 
of two macroscopic atomic ensembles, Physical Review A 83(5): 1-19. 

Nielsen, M. A. & Chuang, I. L. (2011). Quantum Computation and Quan- 
tum Information: 10th Anniversary Edition^ Quantum Computation and 
Quantum Information, Cambridge University Press. 

Oh, S. (2008). Quantum computational method of finding the ground-state 
energy and expectation values. Physical Review A 77(1): 012326. 

Ohzeki, M. & Nishimori, H. (2011). Quantum Annealing: An Introduction and 
New Developments, Journal of Computational and Theoretical Nanoscience 
8(6): 963-971. 

Oliveira, R. & Terhal, B. M. (2005). The complexity of quantum spin systems 
on a two-dimensional square lattice. Quant. Inf, Comp. Vol. 8, No. 10, pp. 
0900-0924 (2008). 

Ortiz, C, Gubernatis, J., Knill, E. & Laflamme, R. (2001). Quantum algorithms 
for fermionic simulations, Physical Review A 64(2): 1-14. 

Osborne, T. J. (2011). Hamiltonian complexity, arXiv:1106.5875 p. 14. 

Papageorgiou, A. & Zhang, C. (2010). On the efficiency of quantum algorithms 
for hamiltonian simulation. Quantum Information Processing pp. 1-21. 

Pastawski, F., Clemente, L. & Cirac, J. (2011). Quantum memories based on 
engineered dissipation. Physical Review A 83(1): 1-12. 

Perdomo, A., Truncik, C, Tubert-Brohman, I., Rose, G. & Aspuru-Guzik, A. 
(2008). Construction of model Hamiltonians for adiabatic quantum com- 
putation and its application to finding low-energy conformations of lattice 
protein models. Physical Review A 78(1): 1-15. 

Pople, J. (1999). Nobel Lecture: Quantum chemical models. Reviews of Modern 
Physics 71(5): 1267-1274. 

Poulin, D. & Wocjan, P. (2009a). Preparing Ground States of Quantum 
Many-Body Systems on a Quantum Computer, Physical Review Letters 
102(13): 1-4. 

Poulin, D. & Wocjan, P. (2009b). Sampling from the Thermal Quantum Gibbs 
State and Evaluating Partition Functions with a Quantum Computer, 
Physical Review Letters 103(22). 

Raesi, S., Wiebe, N. & Sanders, B. C. (2011). Designing Quantum Circuits for 
Efficient Many-Body Quantum Simulation, Arxiv:1108.4318 . 

Rassolov, V. a. & Garashchuk, S. (2008). Computational complexity in quantum 
chemistry. Chemical Physics Letters 464(4-6): 262-264. 



41 



Rebentrost, P., Mohseni, M., Kassal, I., Lloyd, S. & Aspuru-Guzik, A. 
(2009). Environment-assisted quantum transport, New Journal of Physics 
11: 033003. 

Richter, P. (2007). Quantum speedup of classical mixing processes. Physical 
Review A 76(4). 

Riera, A., Gogolin, C. & Eisert, J. (2011). Thermalization in nature and on a 
quantum computer, arXiv:1102.2389 p. 12. 

Runge, E. & Gross, E. K. U. (1984). Density-functional theory for time- 
dependent systems, Phys. Rev. Lett. 52: 997. 

Schuch, N. & Verstraete, F. (2009a). Computational complexity of interacting 
electrons and fundamental limitations of density functional theory, Nature 
Physics 5{10): 732-735. 

Schuch, N. & Verstraete, F. (2009b). Computational Complexity of interacting 
electrons and fundamental limitations of density functional theory. Nature 
Phys. 5: 732. Also see appendix of arxiv:0712.0483. 

Schulman, L., Mor, T. & Weinstein, Y. (2005). Physical Limits of Heat-Bath 
Algorithmic Cooling, Physical Review Letters 94(12): 1-4. 

Shor, P. W. (1997). Polynomial-Time Algorithms for Prime Factorization and 
Discrete Logarithms on a Quantum Computer, SIAM Journal on Comput- 
ing 26(5): 1484. 

Somaroo, S., Tseng, C, Havel, T., Laflamme, R. & Cory, D. (1999). Quan- 
tum Simulations on a Quantum Computer, Physical Review Letters 
82(26): 5381-5384. 

Somma, R., Batista, C. & Ortiz, G. (2007). Quantum Approach to Classical 
Statistical Mechanics, Physical Review Letters 99(3): 1-4. 

Somma, R., Boixo, S., Barnum, H. & Knill, E. (2008). Quantum Simulations of 
Classical Annealing Processes, Physical Review Letters 101(13). 

Somma, R., Ortiz, G., Gubernatis, J. E., Knill, E. & Laflamme, R. (2002). 
Simulating physical phenomena by quantum networks. Physical Review A 
65(4). 

Somma, R., Ortiz, G., Knill, E. & Gubernatis, J. (2003). Quantum Simulations 
of Physics Problems, Proceedings of SPIE 5105: 12. 

Stolze, J. & Sutcr, D. (2008). Quantum computing: a short course from theory 
to experiment, Physics textbook, Wiley- VCH. 

Strini, G. (2002). Error Sensitivity of a Quantum Simulator I: a First Example, 
Fortschritte der Physik 50(2): 171. 



42 



Szegedy, M. (2004). Quantum Speed-Up of Markov Chain Based Algorithms, 
45th Annual IEEE Symposium on Foundations of Computer Science, IEEE, 
pp. 32-41. 

Temme, K., Osborne, T. J., Vohbrecht, K. G., PouUn, D. & Verstraete, F. 
(2011). Quantum Metropohs samphng.. Nature 471(7336): 87-90. 

Tempel, D. G. & Aspuru-Guzik, A. (2011). Quantum computing without wave- 
functions: Time-dependent density functional theory for universal quantum 
computation, arXiv:1108.0097 . 

Terhal, B. M. & DiVincenzo, D. P. (2000). Problem of equilibration and the 
computation of correlation functions on a quantum computer. Physical 
Review A 61(2). 

van Leeuwen, R. (1999). Mapping from Densities to Potentials in Time- 
Dependent Density-Functional Theory, Phys. Rev. Lett. 82(19): 3863-3866. 

Verstraete, F., Wolf, M. M. & Ignacio Cirac, J. (2009). Quantum computa- 
tion and quantum-state engineering driven by dissipation. Nature Physics 
5(9): 633-636. 

Wang, H., Ashhab, S. & Nori, F. (2011). Quantum algorithm for simulating the 
dynamics of an open quantum system. Physical Review A 83(6): 1-11. 

Wang, H., Kais, S., Aspuru-Guzik, A. & Hoffmann, M. R. (2008). Quantum 
algorithm for obtaining the energy spectrum of molecular systems.. Physical 
chemistry chemical physics : PCCP 10(35): 5388-93. 

Ward, N. J., Kassal, I. & Aspuru-Guzik, A. (2009). Preparation of many- 
body states for quantum simulation.. The Journal of chemical physics 
130(19): 194105. 

Waugh, J., Huber, L. & Haeberlen, U. (1968). Approach to high-resolution 
NMR in solids. Physics Review Letters 20(5): 180-182. 

Wei, T.-C., Mosca, M. & Nayak, A. (2010). Interacting Boson Problems Can 
Be QMA Hard, Physical Review Letters 104(4): 1-4. 

Whitfield, J. D., Biamonte, J. & Aspuru-Guzik, A. (2011). Simulation of 
electronic structure Hamiltonians using quantum computers. Molecular 
Physics 109(5): 735-750. 

Wiesner, S. (1996). Simulations of Many-Body Quantum Systems by a Quantum 
Computer, quant-ph/9603028 . 

Williams, C. P. (2010). Explorations in Quantum Computing, Texts in Computer 
Science, Springer. 

Wocjan, P. & Abeyesinghe, A. (2008). Speedup via quantum sampling. Physical 
Review A 78(4). 



43 



Wocjan, P., Chiang, C.-F., Nagaj, D. & Abeyesinghe, A. (2009). Quantum 
algorithm for approximating partition functions, Physical Review A 80(2). 

Wu, L.-a., Byrd, M. & Lidar, D. (2002). Polynomial-Time Simulation of Pairing 
Models on a Quantum Computer, Physical Review Letters 89(5): 1-4. 

Young, a., Knysh, S. & Smclyanskiy, V. (2008). Size Dependence of the Mini- 
mum Excitation Gap in the Quantum Adiabatic Algorithm, Physical Re- 
view Letters 101(17): 1-4. 

Young, a. P. & Smelyanskiy, V. N. (2010). First-Order Phase Transition in the 
Quantum Adiabatic Algorithm, Physical Review Letters 104(2): 1-4. 

Yung, M.-H. & Aspuru-Guzik, A. (2012). A Quantum-Quantum Metropolis 
Algorithm, Proc. Nat. Acad. Sci. In press': 7. 

Yung, M.-H., Boixo, S. & Aspuru-Guzik, A. (2011). Algorithmic quantum cool- 
ing via random walk, under preparation . 

Yung, M.-H., Nagaj, D., Whitfield, J. & Aspuru-Guzik, A. (2010). Simula- 
tion of classical thermal states on a quantum computer: A transfer-matrix 
approach. Physical Review A 82(6): 5. 

Zalka, C. (1998a). Efficient Simulation of Quantum Systems by Quantum Com- 
puters, Fortschritte der Physik 46(6-8): 877-879. 

Zalka, C. (1998b). Simulating quantum systems on a quantum computer. Pro- 
ceedings of the Royal Society A: Mathematical, Physical and Engineering 
Sciences 454(1969): 313-322. 

Zhang, J., Yung, M.-H., Laflamme, R., Aspuru-Guzik, A. & Baugh, J. (2011). 
Digital Quantum Simulation of the Statistical Mechanics of a Frustrated 
Magnet, arXiv:1108.3270 p. 7. 



44 



